ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/rt/ambcomp.c

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.29 by greg, Wed Apr 23 06:04:17 2014 UTC vs.
Revision 2.63 by greg, Thu Jun 19 16:26:55 2014 UTC

#	Line 8 | Line 8 | static const char      RCSid[] = "$Id$";
8		*      for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz
9		*      from ACM SIGGRAPH Asia 2012 conference proceedings.
10		*
11	+	*  Added book-keeping optimization to avoid calculations that would
12	+	*      cancel due to traversal both directions on edges that are adjacent
13	+	*      to same-valued triangles.  This cuts about half of Hessian math.
14	+	*
15		*  Declarations of external symbols in ambient.h
16		*/
17
#	Line 17 | Line 21 | static const char      RCSid[] = "$Id$";
21		#include  "ambient.h"
22		#include  "random.h"
23
24	<	#ifdef NEWAMB
24	>	#ifndef OLDAMB
25
26		extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
27
28		typedef struct {
29	+	COLOR   v;              /* hemisphere sample value */
30	+	float   d;              /* reciprocal distance (1/rt) */
31	+	FVECT   p;              /* intersection point */
32	+	} AMBSAMP;              /* sample value */
33	+
34	+	typedef struct {
35		RAY     *rp;            /* originating ray sample */
26	–	FVECT   ux, uy;         /* tangent axis unit vectors */
36		int     ns;             /* number of samples per axis */
37	+	int     sampOK;         /* acquired full sample set? */
38		COLOR   acoef;          /* division contribution coefficient */
39	<	struct s_ambsamp {
40	<	COLOR   v;              /* hemisphere sample value */
41	<	float   p[3];           /* intersection point */
32	<	} sa[1];                /* sample array (extends struct) */
39	>	double  acol[3];        /* accumulated color */
40	>	FVECT   ux, uy;         /* tangent axis unit vectors */
41	>	AMBSAMP sa[1];          /* sample array (extends struct) */
42		}  AMBHEMI;             /* ambient sample hemisphere */
43
44	<	#define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
44	>	#define AI(h,i,j)       ((i)*(h)->ns + (j))
45	>	#define ambsam(h,i,j)   (h)->sa[AI(h,i,j)]
46
47		typedef struct {
48	<	FVECT   r_i, r_i1, e_i;
49	<	double  nf, I1, I2, J2;
48	>	FVECT   r_i, r_i1, e_i, rcp, rI2_eJ2;
49	>	double  I1, I2;
50		} FFTRI;                /* vectors and coefficients for Hessian calculation */
51
52
53	+	static int
54	+	ambsample(                              /* initial ambient division sample */
55	+	AMBHEMI *hp,
56	+	int     i,
57	+	int     j,
58	+	int     n
59	+	)
60	+	{
61	+	AMBSAMP *ap = &ambsam(hp,i,j);
62	+	RAY     ar;
63	+	int     hlist[3], ii;
64	+	double  spt[2], zd;
65	+	/* generate hemispherical sample */
66	+	/* ambient coefficient for weight */
67	+	if (ambacc > FTINY)
68	+	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
69	+	else
70	+	copycolor(ar.rcoef, hp->acoef);
71	+	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
72	+	return(0);
73	+	if (ambacc > FTINY) {
74	+	multcolor(ar.rcoef, hp->acoef);
75	+	scalecolor(ar.rcoef, 1./AVGREFL);
76	+	}
77	+	hlist[0] = hp->rp->rno;
78	+	hlist[1] = j;
79	+	hlist[2] = i;
80	+	multisamp(spt, 2, urand(ilhash(hlist,3)+n));
81	+	if (!n) {                       /* avoid border samples for n==0 */
82	+	if ((spt[0] < 0.1) | (spt[0] >= 0.9))
83	+	spt[0] = 0.1 + 0.8*frandom();
84	+	if ((spt[1] < 0.1) | (spt[1] >= 0.9))
85	+	spt[1] = 0.1 + 0.8*frandom();
86	+	}
87	+	SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
88	+	zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
89	+	for (ii = 3; ii--; )
90	+	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
91	+	spt[1]*hp->uy[ii] +
92	+	zd*hp->rp->ron[ii];
93	+	checknorm(ar.rdir);
94	+	dimlist[ndims++] = AI(hp,i,j) + 90171;
95	+	rayvalue(&ar);                  /* evaluate ray */
96	+	ndims--;
97	+	if (ar.rt <= FTINY)
98	+	return(0);              /* should never happen */
99	+	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
100	+	if (ar.rt*ap->d < 1.0)          /* new/closer distance? */
101	+	ap->d = 1.0/ar.rt;
102	+	if (!n) {                       /* record first vertex & value */
103	+	if (ar.rt > 10.0*thescene.cusize)
104	+	ar.rt = 10.0*thescene.cusize;
105	+	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
106	+	copycolor(ap->v, ar.rcol);
107	+	} else {                        /* else update recorded value */
108	+	hp->acol[RED] -= colval(ap->v,RED);
109	+	hp->acol[GRN] -= colval(ap->v,GRN);
110	+	hp->acol[BLU] -= colval(ap->v,BLU);
111	+	zd = 1.0/(double)(n+1);
112	+	scalecolor(ar.rcol, zd);
113	+	zd *= (double)n;
114	+	scalecolor(ap->v, zd);
115	+	addcolor(ap->v, ar.rcol);
116	+	}
117	+	addcolor(hp->acol, ap->v);      /* add to our sum */
118	+	return(1);
119	+	}
120	+
121	+
122	+	/* Estimate errors based on ambient division differences */
123	+	static float *
124	+	getambdiffs(AMBHEMI *hp)
125	+	{
126	+	float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
127	+	float   *ep;
128	+	AMBSAMP *ap;
129	+	double  b, d2;
130	+	int     i, j;
131	+
132	+	if (earr == NULL)               /* out of memory? */
133	+	return(NULL);
134	+	/* compute squared neighbor diffs */
135	+	for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
136	+	for (j = 0; j < hp->ns; j++, ap++, ep++) {
137	+	b = bright(ap[0].v);
138	+	if (i) {                /* from above */
139	+	d2 = b - bright(ap[-hp->ns].v);
140	+	d2 *= d2;
141	+	ep[0] += d2;
142	+	ep[-hp->ns] += d2;
143	+	}
144	+	if (!j) continue;
145	+	/* from behind */
146	+	d2 = b - bright(ap[-1].v);
147	+	d2 *= d2;
148	+	ep[0] += d2;
149	+	ep[-1] += d2;
150	+	if (!i) continue;
151	+	/* diagonal */
152	+	d2 = b - bright(ap[-hp->ns-1].v);
153	+	d2 *= d2;
154	+	ep[0] += d2;
155	+	ep[-hp->ns-1] += d2;
156	+	}
157	+	/* correct for number of neighbors */
158	+	earr[0] *= 8./3.;
159	+	earr[hp->ns-1] *= 8./3.;
160	+	earr[(hp->ns-1)*hp->ns] *= 8./3.;
161	+	earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 8./3.;
162	+	for (i = 1; i < hp->ns-1; i++) {
163	+	earr[i*hp->ns] *= 8./5.;
164	+	earr[i*hp->ns + hp->ns-1] *= 8./5.;
165	+	}
166	+	for (j = 1; j < hp->ns-1; j++) {
167	+	earr[j] *= 8./5.;
168	+	earr[(hp->ns-1)*hp->ns + j] *= 8./5.;
169	+	}
170	+	return(earr);
171	+	}
172	+
173	+
174	+	/* Perform super-sampling on hemisphere (introduces bias) */
175	+	static void
176	+	ambsupersamp(AMBHEMI *hp, int cnt)
177	+	{
178	+	float   *earr = getambdiffs(hp);
179	+	double  e2rem = 0;
180	+	AMBSAMP *ap;
181	+	RAY     ar;
182	+	float   *ep;
183	+	int     i, j, n, nss;
184	+
185	+	if (earr == NULL)               /* just skip calc. if no memory */
186	+	return;
187	+	/* accumulate estimated variances */
188	+	for (ep = earr + hp->ns*hp->ns; ep > earr; )
189	+	e2rem += *--ep;
190	+	ep = earr;                      /* perform super-sampling */
191	+	for (ap = hp->sa, i = 0; i < hp->ns; i++)
192	+	for (j = 0; j < hp->ns; j++, ap++) {
193	+	if (e2rem <= FTINY)
194	+	goto done;      /* nothing left to do */
195	+	nss = *ep/e2rem*cnt + frandom();
196	+	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
197	+	--cnt;
198	+	e2rem -= *ep++;         /* update remainder */
199	+	}
200	+	done:
201	+	free(earr);
202	+	}
203	+
204	+
205		static AMBHEMI *
206	<	inithemi(                       /* initialize sampling hemisphere */
207	<	COLOR   ac,
206	>	samp_hemi(                              /* sample indirect hemisphere */
207	>	COLOR   rcol,
208		RAY     *r,
209		double  wt
210		)
211		{
212		AMBHEMI *hp;
213		double  d;
214	<	int     n, i;
214	>	int     n, i, j;
215		/* set number of divisions */
216		if (ambacc <= FTINY &&
217	<	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
217	>	wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight)))
218		wt = d;                 /* avoid ray termination */
219		n = sqrt(ambdiv * wt) + 0.5;
220		i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
221		if (n < i)
222		n = i;
223		/* allocate sampling array */
224	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
63	<	sizeof(struct s_ambsamp)*(n*n - 1));
224	>	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
225		if (hp == NULL)
226	<	return(NULL);
226	>	error(SYSTEM, "out of memory in samp_hemi");
227		hp->rp = r;
228		hp->ns = n;
229	+	hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0;
230	+	memset(hp->sa, 0, sizeof(AMBSAMP)*n*n);
231	+	hp->sampOK = 0;
232		/* assign coefficient */
233	<	copycolor(hp->acoef, ac);
233	>	copycolor(hp->acoef, rcol);
234		d = 1.0/(n*n);
235		scalecolor(hp->acoef, d);
236		/* make tangent plane axes */
237	<	hp->uy[0] = 0.1 - 0.2*frandom();
238	<	hp->uy[1] = 0.1 - 0.2*frandom();
239	<	hp->uy[2] = 0.1 - 0.2*frandom();
240	<	for (i = 0; i < 3; i++)
241	<	if (r->ron[i] < 0.6 && r->ron[i] > -0.6)
237	>	hp->uy[0] = 0.5 - frandom();
238	>	hp->uy[1] = 0.5 - frandom();
239	>	hp->uy[2] = 0.5 - frandom();
240	>	for (i = 3; i--; )
241	>	if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
242		break;
243	<	if (i >= 3)
244	<	error(CONSISTENCY, "bad ray direction in inithemi()");
243	>	if (i < 0)
244	>	error(CONSISTENCY, "bad ray direction in samp_hemi");
245		hp->uy[i] = 1.0;
246		VCROSS(hp->ux, hp->uy, r->ron);
247		normalize(hp->ux);
248		VCROSS(hp->uy, r->ron, hp->ux);
249	<	/* we're ready to sample */
250	<	return(hp);
249	>	/* sample divisions */
250	>	for (i = hp->ns; i--; )
251	>	for (j = hp->ns; j--; )
252	>	hp->sampOK += ambsample(hp, i, j, 0);
253	>	copycolor(rcol, hp->acol);
254	>	if (!hp->sampOK) {              /* utter failure? */
255	>	free(hp);
256	>	return(NULL);
257	>	}
258	>	if (hp->sampOK < hp->ns*hp->ns) {
259	>	hp->sampOK *= -1;       /* soft failure */
260	>	return(hp);
261	>	}
262	>	n = ambssamp*wt + 0.5;
263	>	if (n > 8) {                    /* perform super-sampling? */
264	>	ambsupersamp(hp, n);
265	>	copycolor(rcol, hp->acol);
266	>	}
267	>	return(hp);                     /* all is well */
268		}
269
270
271	<	static struct s_ambsamp *
272	<	ambsample(                              /* sample an ambient direction */
273	<	AMBHEMI *hp,
93	<	int     i,
94	<	int     j
95	<	)
271	>	/* Return brightness of farthest ambient sample */
272	>	static double
273	>	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
274		{
275	<	struct s_ambsamp        *ap = &ambsamp(hp,i,j);
276	<	RAY                     ar;
277	<	double                  spt[2], zd;
278	<	int                     ii;
101	<	/* ambient coefficient for weight */
102	<	if (ambacc > FTINY)
103	<	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
104	<	else
105	<	copycolor(ar.rcoef, hp->acoef);
106	<	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) {
107	<	setcolor(ap->v, 0., 0., 0.);
108	<	VCOPY(ap->p, hp->rp->rop);
109	<	return(NULL);           /* no sample taken */
275	>	if (hp->sa[n1].d <= hp->sa[n2].d) {
276	>	if (hp->sa[n1].d <= hp->sa[n3].d)
277	>	return(colval(hp->sa[n1].v,CIEY));
278	>	return(colval(hp->sa[n3].v,CIEY));
279		}
280	<	if (ambacc > FTINY) {
281	<	multcolor(ar.rcoef, hp->acoef);
282	<	scalecolor(ar.rcoef, 1./AVGREFL);
114	<	}
115	<	/* generate hemispherical sample */
116	<	SDsquare2disk(spt,      (i+.1+.8*frandom())/hp->ns,
117	<	(j+.1+.8*frandom())/hp->ns );
118	<	zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
119	<	for (ii = 3; ii--; )
120	<	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
121	<	spt[1]*hp->uy[ii] +
122	<	zd*hp->rp->ron[ii];
123	<	checknorm(ar.rdir);
124	<	dimlist[ndims++] = i*hp->ns + j + 90171;
125	<	rayvalue(&ar);                  /* evaluate ray */
126	<	ndims--;
127	<	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
128	<	copycolor(ap->v, ar.rcol);
129	<	if (ar.rt > 20.0*maxarad)       /* limit vertex distance */
130	<	VSUM(ap->p, ar.rorg, ar.rdir, 20.0*maxarad);
131	<	else
132	<	VCOPY(ap->p, ar.rop);
133	<	return(ap);
280	>	if (hp->sa[n2].d <= hp->sa[n3].d)
281	>	return(colval(hp->sa[n2].v,CIEY));
282	>	return(colval(hp->sa[n3].v,CIEY));
283		}
284
285
286		/* Compute vectors and coefficients for Hessian/gradient calcs */
287		static void
288	<	comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop)
288	>	comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1)
289		{
290	<	FVECT   v1;
291	<	double  dot_e, dot_er, dot_r, dot_r1;
290	>	double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
291	>	int     ii;
292
293	<	VSUB(ftp->r_i, ap0, rop);
294	<	VSUB(ftp->r_i1, ap1, rop);
295	<	VSUB(ftp->e_i, ap1, ap0);
296	<	VCROSS(v1, ftp->e_i, ftp->r_i);
297	<	ftp->nf = 1.0/DOT(v1,v1);
293	>	VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop);
294	>	VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop);
295	>	VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p);
296	>	VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
297	>	rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
298		dot_e = DOT(ftp->e_i,ftp->e_i);
299		dot_er = DOT(ftp->e_i, ftp->r_i);
300	<	dot_r = DOT(ftp->r_i,ftp->r_i);
301	<	dot_r1 = DOT(ftp->r_i1,ftp->r_i1);
302	<	ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) / sqrt(dot_r*dot_r1) ) *
303	<	sqrt( ftp->nf );
304	<	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r +
305	<	dot_e*ftp->I1 )*0.5*ftp->nf;
306	<	ftp->J2 =  0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) -
307	<	dot_er/dot_e*ftp->I2;
300	>	rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i);
301	>	rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
302	>	ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) *
303	>	sqrt( rdot_cp );
304	>	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
305	>	dot_e*ftp->I1 )*0.5*rdot_cp;
306	>	J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
307	>	for (ii = 3; ii--; )
308	>	ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
309		}
310
311
#	Line 176 | Line 326 | compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
326		static void
327		comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
328		{
329	<	FVECT   v1, v2;
329	>	FVECT   ncp;
330		FVECT   m1[3], m2[3], m3[3], m4[3];
331		double  d1, d2, d3, d4;
332		double  I3, J3, K3;
#	Line 186 | Line 336 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
336		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
337		d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
338		d4 = DOT(ftp->e_i, ftp->r_i);
339	<	I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 +
340	<	3.0/d3*ftp->I2 );
339	>	I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 )
340	>	/ ( 4.0*DOT(ftp->rcp,ftp->rcp) );
341		J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3;
342		K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3);
343		/* intermediate matrices */
344	<	VCROSS(v1, nrm, ftp->e_i);
345	<	for (j = 3; j--; )
196	<	v2[j] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j];
197	<	compose_matrix(m1, v1, v2);
344	>	VCROSS(ncp, nrm, ftp->e_i);
345	>	compose_matrix(m1, ncp, ftp->rI2_eJ2);
346		compose_matrix(m2, ftp->r_i, ftp->r_i);
347		compose_matrix(m3, ftp->e_i, ftp->e_i);
348		compose_matrix(m4, ftp->r_i, ftp->e_i);
349	<	VCROSS(v1, ftp->r_i, ftp->e_i);
202	<	d1 = DOT(nrm, v1);
349	>	d1 = DOT(nrm, ftp->rcp);
350		d2 = -d1*ftp->I2;
351		d1 *= 2.0;
352		for (i = 3; i--; )              /* final matrix sum */
#	Line 229 | Line 376 | rev_hessian(FVECT hess[3])
376		/* Add to radiometric Hessian from the given triangle */
377		static void
378		add2hessian(FVECT hess[3], FVECT ehess1[3],
379	<	FVECT ehess2[3], FVECT ehess3[3], COLORV v)
379	>	FVECT ehess2[3], FVECT ehess3[3], double v)
380		{
381		int     i, j;
382
#	Line 243 | Line 390 | add2hessian(FVECT hess[3], FVECT ehess1[3],
390		static void
391		comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
392		{
393	<	FVECT   vcp;
393	>	FVECT   ncp;
394		double  f1;
395		int     i;
396
397	<	VCROSS(vcp, ftp->r_i, ftp->r_i1);
398	<	f1 = 2.0*DOT(nrm, vcp);
252	<	VCROSS(vcp, nrm, ftp->e_i);
397	>	f1 = 2.0*DOT(nrm, ftp->rcp);
398	>	VCROSS(ncp, nrm, ftp->e_i);
399		for (i = 3; i--; )
400	<	grad[i] = (0.5/PI)*( ftp->I1*vcp[i] +
255	<	f1*(ftp->I2*ftp->r_i[i] + ftp->J2*ftp->e_i[i]) );
400	>	grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
401		}
402
403
#	Line 268 | Line 413 | rev_gradient(FVECT grad)
413
414		/* Add to displacement gradient from the given triangle */
415		static void
416	<	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
416	>	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
417		{
418		int     i;
419
#	Line 277 | Line 422 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
422		}
423
424
280	–	/* Return brightness of furthest ambient sample */
281	–	static COLORV
282	–	back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2,
283	–	struct s_ambsamp *ap3, FVECT orig)
284	–	{
285	–	COLORV  vback;
286	–	FVECT   vec;
287	–	double  d2, d2best;
288	–
289	–	VSUB(vec, ap1->p, orig);
290	–	d2best = DOT(vec,vec);
291	–	vback = colval(ap1->v,CIEY);
292	–	VSUB(vec, ap2->p, orig);
293	–	d2 = DOT(vec,vec);
294	–	if (d2 > d2best) {
295	–	d2best = d2;
296	–	vback = colval(ap2->v,CIEY);
297	–	}
298	–	VSUB(vec, ap3->p, orig);
299	–	d2 = DOT(vec,vec);
300	–	if (d2 > d2best)
301	–	return(colval(ap3->v,CIEY));
302	–	return(vback);
303	–	}
304	–
305	–
425		/* Compute anisotropic radii and eigenvector directions */
426	<	static int
426	>	static void
427		eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
428		{
429		double  hess2[2][2];
#	Line 320 | Line 439 | eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
439		hess2[0][1] = DOT(uv[0], b);
440		hess2[1][0] = DOT(uv[1], a);
441		hess2[1][1] = DOT(uv[1], b);
442	<	/* compute eigenvalues */
443	<	if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
444	<	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 ||
445	<	(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY ||
446	<	(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY )
447	<	error(INTERNAL, "bad eigenvalue calculation");
448	<
442	>	/* compute eigenvalue(s) */
443	>	i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
444	>	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]);
445	>	if (i == 1)                     /* double-root (circle) */
446	>	evalue[1] = evalue[0];
447	>	if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) |
448	>	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) {
449	>	ra[0] = ra[1] = maxarad;
450	>	return;
451	>	}
452		if (evalue[0] > evalue[1]) {
453		ra[0] = sqrt(sqrt(4.0/evalue[0]));
454		ra[1] = sqrt(sqrt(4.0/evalue[1]));
#	Line 384 | Line 506 | ambHessian(                            /* anisotropic radii & pos. gradient */
506		}
507		/* compute first row of edges */
508		for (j = 0; j < hp->ns-1; j++) {
509	<	comp_fftri(&fftr, ambsamp(hp,0,j).p,
388	<	ambsamp(hp,0,j+1).p, hp->rp->rop);
509	>	comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1));
510		if (hessrow != NULL)
511		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
512		if (gradrow != NULL)
#	Line 395 | Line 516 | ambHessian(                            /* anisotropic radii & pos. gradient */
516		for (i = 0; i < hp->ns-1; i++) {
517		FVECT       hesscol[3];     /* compute first vertical edge */
518		FVECT       gradcol;
519	<	comp_fftri(&fftr, ambsamp(hp,i,0).p,
399	<	ambsamp(hp,i+1,0).p, hp->rp->rop);
519	>	comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0));
520		if (hessrow != NULL)
521		comp_hessian(hesscol, &fftr, hp->rp->ron);
522		if (gradrow != NULL)
#	Line 404 | Line 524 | ambHessian(                            /* anisotropic radii & pos. gradient */
524		for (j = 0; j < hp->ns-1; j++) {
525		FVECT   hessdia[3];     /* compute triangle contributions */
526		FVECT   graddia;
527	<	COLORV  backg;
528	<	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
529	<	&ambsamp(hp,i+1,j), hp->rp->rop);
527	>	double  backg;
528	>	backg = back_ambval(hp, AI(hp,i,j),
529	>	AI(hp,i,j+1), AI(hp,i+1,j));
530		/* diagonal (inner) edge */
531	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
412	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
531	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j));
532		if (hessrow != NULL) {
533		comp_hessian(hessdia, &fftr, hp->rp->ron);
534		rev_hessian(hesscol);
535		add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
536		}
537	<	if (gradient != NULL) {
537	>	if (gradrow != NULL) {
538		comp_gradient(graddia, &fftr, hp->rp->ron);
539		rev_gradient(gradcol);
540		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
541		}
542		/* initialize edge in next row */
543	<	comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
425	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
543	>	comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j));
544		if (hessrow != NULL)
545		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
546		if (gradrow != NULL)
547		comp_gradient(gradrow[j], &fftr, hp->rp->ron);
548		/* new column edge & paired triangle */
549	<	backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
550	<	&ambsamp(hp,i+1,j), hp->rp->rop);
551	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
434	<	hp->rp->rop);
549	>	backg = back_ambval(hp, AI(hp,i+1,j+1),
550	>	AI(hp,i+1,j), AI(hp,i,j+1));
551	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1));
552		if (hessrow != NULL) {
553		comp_hessian(hesscol, &fftr, hp->rp->ron);
554		rev_hessian(hessdia);
#	Line 454 | Line 571 | ambHessian(                            /* anisotropic radii & pos. gradient */
571
572		if (ra != NULL)                 /* extract eigenvectors & radii */
573		eigenvectors(uv, ra, hessian);
574	<	if (pg != NULL) {               /* tangential position gradient/PI */
575	<	pg[0] = DOT(gradient, uv[0]) / PI;
576	<	pg[1] = DOT(gradient, uv[1]) / PI;
574	>	if (pg != NULL) {               /* tangential position gradient */
575	>	pg[0] = DOT(gradient, uv[0]);
576	>	pg[1] = DOT(gradient, uv[1]);
577		}
578		}
579
#	Line 465 | Line 582 | ambHessian(                            /* anisotropic radii & pos. gradient */
582		static void
583		ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
584		{
585	<	struct s_ambsamp        *ap;
586	<	double                  dgsum[2];
587	<	int                     n;
588	<	FVECT                   vd;
589	<	double                  gfact;
585	>	AMBSAMP *ap;
586	>	double  dgsum[2];
587	>	int     n;
588	>	FVECT   vd;
589	>	double  gfact;
590
591		dgsum[0] = dgsum[1] = 0.0;      /* sum values times -tan(theta) */
592		for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
#	Line 477 | Line 594 | ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
594		VSUB(vd, ap->p, hp->rp->rop);
595		/* brightness over cosine factor */
596		gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
597	<	/* -sine = -proj_radius/vd_length */
598	<	dgsum[0] += DOT(uv[1], vd) * gfact;
599	<	dgsum[1] -= DOT(uv[0], vd) * gfact;
597	>	/* sine = proj_radius/vd_length */
598	>	dgsum[0] -= DOT(uv[1], vd) * gfact;
599	>	dgsum[1] += DOT(uv[0], vd) * gfact;
600		}
601		dg[0] = dgsum[0] / (hp->ns*hp->ns);
602		dg[1] = dgsum[1] / (hp->ns*hp->ns);
603		}
604
605
606	+	/* Compute potential light leak direction flags for cache value */
607	+	static uint32
608	+	ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
609	+	{
610	+	const double    max_d = 1.0/(minarad*ambacc + 0.001);
611	+	const double    ang_res = 0.5*PI/(hp->ns-1);
612	+	const double    ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY));
613	+	double          avg_d = 0;
614	+	uint32          flgs = 0;
615	+	FVECT           vec;
616	+	double          u, v;
617	+	double          ang, a1;
618	+	int             i, j;
619	+	/* don't bother for a few samples */
620	+	if (hp->ns < 12)
621	+	return(0);
622	+	/* check distances overhead */
623	+	for (i = hp->ns*3/4; i-- > hp->ns>>2; )
624	+	for (j = hp->ns*3/4; j-- > hp->ns>>2; )
625	+	avg_d += ambsam(hp,i,j).d;
626	+	avg_d *= 4.0/(hp->ns*hp->ns);
627	+	if (avg_d*r0 >= 1.0)            /* ceiling too low for corral? */
628	+	return(0);
629	+	if (avg_d >= max_d)             /* insurance */
630	+	return(0);
631	+	/* else circle around perimeter */
632	+	for (i = 0; i < hp->ns; i++)
633	+	for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) {
634	+	AMBSAMP *ap = &ambsam(hp,i,j);
635	+	if ((ap->d <= FTINY) | (ap->d >= max_d))
636	+	continue;       /* too far or too near */
637	+	VSUB(vec, ap->p, hp->rp->rop);
638	+	u = DOT(vec, uv[0]);
639	+	v = DOT(vec, uv[1]);
640	+	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v)
641	+	continue;       /* occluder outside ellipse */
642	+	ang = atan2a(v, u);     /* else set direction flags */
643	+	for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step)
644	+	flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0)));
645	+	}
646	+	/* add low-angle incident (< 20deg) */
647	+	if (fabs(hp->rp->rod) <= 0.342) {
648	+	u = -DOT(hp->rp->rdir, uv[0]);
649	+	v = -DOT(hp->rp->rdir, uv[1]);
650	+	if ((r0*r0*u*u + r1*r1*v*v) > hp->rp->rot*hp->rp->rot) {
651	+	ang = atan2a(v, u);
652	+	ang += 2.*PI*(ang < 0);
653	+	ang *= 16/PI;
654	+	if ((ang < .5) | (ang >= 31.5))
655	+	flgs |= 0x80000001;
656	+	else
657	+	flgs |= 3L<<(int)(ang-.5);
658	+	}
659	+	}
660	+	return(flgs);
661	+	}
662	+
663	+
664		int
665		doambient(                              /* compute ambient component */
666		COLOR   rcol,                   /* input/output color */
#	Line 494 | Line 669 | doambient(                             /* compute ambient component */
669		FVECT   uv[2],                  /* returned (optional) */
670		float   ra[2],                  /* returned (optional) */
671		float   pg[2],                  /* returned (optional) */
672	<	float   dg[2]                   /* returned (optional) */
672	>	float   dg[2],                  /* returned (optional) */
673	>	uint32  *crlp                   /* returned (optional) */
674		)
675		{
676	<	AMBHEMI                 *hp = inithemi(rcol, r, wt);
677	<	int                     cnt = 0;
678	<	FVECT                   my_uv[2];
679	<	double                  d, acol[3];
680	<	struct s_ambsamp        *ap;
681	<	int                     i, j;
506	<	/* check/initialize */
507	<	if (hp == NULL)
508	<	return(0);
676	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
677	>	FVECT   my_uv[2];
678	>	double  d, K;
679	>	AMBSAMP *ap;
680	>	int     i;
681	>	/* clear return values */
682		if (uv != NULL)
683		memset(uv, 0, sizeof(FVECT)*2);
684		if (ra != NULL)
#	Line 514 | Line 687 | doambient(                             /* compute ambient component */
687		pg[0] = pg[1] = 0.0;
688		if (dg != NULL)
689		dg[0] = dg[1] = 0.0;
690	<	/* sample the hemisphere */
691	<	acol[0] = acol[1] = acol[2] = 0.0;
692	<	for (i = hp->ns; i--; )
693	<	for (j = hp->ns; j--; )
694	<	if ((ap = ambsample(hp, i, j)) != NULL) {
695	<	addcolor(acol, ap->v);
696	<	++cnt;
697	<	}
698	<	if (!cnt) {
526	<	setcolor(rcol, 0.0, 0.0, 0.0);
527	<	free(hp);
528	<	return(0);              /* no valid samples */
690	>	if (crlp != NULL)
691	>	*crlp = 0;
692	>	if (hp == NULL)                 /* sampling falure? */
693	>	return(0);
694	>
695	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) ||
696	>	(hp->sampOK < 0) | (hp->ns < 4)) {
697	>	free(hp);               /* Hessian not requested/possible */
698	>	return(-1);             /* value-only return value */
699		}
700	<	copycolor(rcol, acol);          /* final indirect irradiance/PI */
701	<	if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
702	<	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
703	<	free(hp);
704	<	return(-1);             /* no radius or gradient calc. */
700	>	if ((d = bright(rcol)) > FTINY) {       /* normalize Y values */
701	>	d = 0.99*(hp->ns*hp->ns)/d;
702	>	K = 0.01;
703	>	} else {                        /* or fall back on geometric Hessian */
704	>	K = 1.0;
705	>	pg = NULL;
706	>	dg = NULL;
707	>	crlp = NULL;
708		}
536	–	multcolor(acol, hp->acoef);     /* normalize Y values */
537	–	if ((d = bright(acol)) > FTINY)
538	–	d = 1.0/d;
539	–	else
540	–	d = 0.0;
709		ap = hp->sa;                    /* relative Y channel from here on... */
710		for (i = hp->ns*hp->ns; i--; ap++)
711	<	colval(ap->v,CIEY) = bright(ap->v)*d + 0.0314;
711	>	colval(ap->v,CIEY) = bright(ap->v)*d + K;
712
713		if (uv == NULL)                 /* make sure we have axis pointers */
714		uv = my_uv;
#	Line 551 | Line 719 | doambient(                             /* compute ambient component */
719		ambdirgrad(hp, uv, dg);
720
721		if (ra != NULL) {               /* scale/clamp radii */
722	+	if (pg != NULL) {
723	+	if (ra[0]*(d = fabs(pg[0])) > 1.0)
724	+	ra[0] = 1.0/d;
725	+	if (ra[1]*(d = fabs(pg[1])) > 1.0)
726	+	ra[1] = 1.0/d;
727	+	if (ra[0] > ra[1])
728	+	ra[0] = ra[1];
729	+	}
730		if (ra[0] < minarad) {
731		ra[0] = minarad;
732		if (ra[1] < minarad)
733		ra[1] = minarad;
734		}
735	<	ra[0] *= d = 1.0/sqrt(sqrt(wt));
735	>	ra[0] *= d = 1.0/sqrt(wt);
736		if ((ra[1] *= d) > 2.0*ra[0])
737		ra[1] = 2.0*ra[0];
738		if (ra[1] > maxarad) {
739		ra[1] = maxarad;
740		if (ra[0] > maxarad)
741		ra[0] = maxarad;
742	+	}
743	+	/* flag encroached directions */
744	+	if ((wt >= 0.89*AVGREFL) & (crlp != NULL))
745	+	*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc);
746	+	if (pg != NULL) {       /* cap gradient if necessary */
747	+	d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];
748	+	if (d > 1.0) {
749	+	d = 1.0/sqrt(d);
750	+	pg[0] *= d;
751	+	pg[1] *= d;
752	+	}
753		}
754		}
755		free(hp);                       /* clean up and return */

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)