[ViewVC] Diff of: radiance/ray/src/rt/ambcomp.c

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.27 by greg, Sat Apr 19 02:39:44 2014 UTC vs.
Revision 2.64 by greg, Tue Aug 19 15:04:40 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	+	/* avoid coincident samples */
82	+	if (!n && (0 < i) & (i < hp->ns-1) &&
83	+	(0 < j) & (j < hp->ns-1)) {
84	+	if ((spt[0] < 0.1) \| (spt[0] >= 0.9))
85	+	spt[0] = 0.1 + 0.8*frandom();
86	+	if ((spt[1] < 0.1) \| (spt[1] >= 0.9))
87	+	spt[1] = 0.1 + 0.8*frandom();
88	+	}
89	+	SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
90	+	zd = sqrt(1. - spt[0]spt[0] - spt[1]spt[1]);
91	+	for (ii = 3; ii--; )
92	+	ar.rdir[ii] = spt[0]*hp->ux[ii] +
93	+	spt[1]*hp->uy[ii] +
94	+	zd*hp->rp->ron[ii];
95	+	checknorm(ar.rdir);
96	+	dimlist[ndims++] = AI(hp,i,j) + 90171;
97	+	rayvalue(&ar); /* evaluate ray */
98	+	ndims--;
99	+	if (ar.rt <= FTINY)
100	+	return(0); /* should never happen */
101	+	multcolor(ar.rcol, ar.rcoef); /* apply coefficient */
102	+	if (ar.rtap->d < 1.0) / new/closer distance? */
103	+	ap->d = 1.0/ar.rt;
104	+	if (!n) { /* record first vertex & value */
105	+	if (ar.rt > 10.0*thescene.cusize)
106	+	ar.rt = 10.0*thescene.cusize;
107	+	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
108	+	copycolor(ap->v, ar.rcol);
109	+	} else { /* else update recorded value */
110	+	hp->acol[RED] -= colval(ap->v,RED);
111	+	hp->acol[GRN] -= colval(ap->v,GRN);
112	+	hp->acol[BLU] -= colval(ap->v,BLU);
113	+	zd = 1.0/(double)(n+1);
114	+	scalecolor(ar.rcol, zd);
115	+	zd *= (double)n;
116	+	scalecolor(ap->v, zd);
117	+	addcolor(ap->v, ar.rcol);
118	+	}
119	+	addcolor(hp->acol, ap->v); /* add to our sum */
120	+	return(1);
121	+	}
122	+
123	+
124	+	/* Estimate errors based on ambient division differences */
125	+	static float *
126	+	getambdiffs(AMBHEMI *hp)
127	+	{
128	+	float earr = (float )calloc(hp->ns*hp->ns, sizeof(float));
129	+	float *ep;
130	+	AMBSAMP *ap;
131	+	double b, d2;
132	+	int i, j;
133	+
134	+	if (earr == NULL) /* out of memory? */
135	+	return(NULL);
136	+	/* compute squared neighbor diffs */
137	+	for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
138	+	for (j = 0; j < hp->ns; j++, ap++, ep++) {
139	+	b = bright(ap[0].v);
140	+	if (i) { /* from above */
141	+	d2 = b - bright(ap[-hp->ns].v);
142	+	d2 *= d2;
143	+	ep[0] += d2;
144	+	ep[-hp->ns] += d2;
145	+	}
146	+	if (!j) continue;
147	+	/* from behind */
148	+	d2 = b - bright(ap[-1].v);
149	+	d2 *= d2;
150	+	ep[0] += d2;
151	+	ep[-1] += d2;
152	+	if (!i) continue;
153	+	/* diagonal */
154	+	d2 = b - bright(ap[-hp->ns-1].v);
155	+	d2 *= d2;
156	+	ep[0] += d2;
157	+	ep[-hp->ns-1] += d2;
158	+	}
159	+	/* correct for number of neighbors */
160	+	earr[0] *= 8./3.;
161	+	earr[hp->ns-1] *= 8./3.;
162	+	earr[(hp->ns-1)hp->ns] = 8./3.;
163	+	earr[(hp->ns-1)hp->ns + hp->ns-1] = 8./3.;
164	+	for (i = 1; i < hp->ns-1; i++) {
165	+	earr[ihp->ns] = 8./5.;
166	+	earr[ihp->ns + hp->ns-1] = 8./5.;
167	+	}
168	+	for (j = 1; j < hp->ns-1; j++) {
169	+	earr[j] *= 8./5.;
170	+	earr[(hp->ns-1)hp->ns + j] = 8./5.;
171	+	}
172	+	return(earr);
173	+	}
174	+
175	+
176	+	/* Perform super-sampling on hemisphere (introduces bias) */
177	+	static void
178	+	ambsupersamp(AMBHEMI *hp, int cnt)
179	+	{
180	+	float *earr = getambdiffs(hp);
181	+	double e2rem = 0;
182	+	AMBSAMP *ap;
183	+	RAY ar;
184	+	float *ep;
185	+	int i, j, n, nss;
186	+
187	+	if (earr == NULL) /* just skip calc. if no memory */
188	+	return;
189	+	/* accumulate estimated variances */
190	+	for (ep = earr + hp->ns*hp->ns; ep > earr; )
191	+	e2rem += *--ep;
192	+	ep = earr; /* perform super-sampling */
193	+	for (ap = hp->sa, i = 0; i < hp->ns; i++)
194	+	for (j = 0; j < hp->ns; j++, ap++) {
195	+	if (e2rem <= FTINY)
196	+	goto done; /* nothing left to do */
197	+	nss = ep/e2remcnt + frandom();
198	+	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
199	+	--cnt;
200	+	e2rem -= ep++; / update remainder */
201	+	}
202	+	done:
203	+	free(earr);
204	+	}
205	+
206	+
207		static AMBHEMI *
208	<	inithemi( /* initialize sampling hemisphere */
209	<	COLOR ac,
208	>	samp_hemi( /* sample indirect hemisphere */
209	>	COLOR rcol,
210		RAY *r,
211		double wt
212		)
213		{
214		AMBHEMI *hp;
215		double d;
216	<	int n, i;
216	>	int n, i, j;
217		/* set number of divisions */
218		if (ambacc <= FTINY &&
219	<	wt > (d = 0.8intens(ac)r->rweight/(ambdiv*minweight)))
219	>	wt > (d = 0.8intens(rcol)r->rweight/(ambdiv*minweight)))
220		wt = d; /* avoid ray termination */
221		n = sqrt(ambdiv * wt) + 0.5;
222		i = 1 + 5(ambacc > FTINY); / minimum number of samples */
223		if (n < i)
224		n = i;
225		/* allocate sampling array */
226	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
63	<	sizeof(struct s_ambsamp)(nn - 1));
226	>	hp = (AMBHEMI )malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)(n*n - 1));
227		if (hp == NULL)
228	<	return(NULL);
228	>	error(SYSTEM, "out of memory in samp_hemi");
229		hp->rp = r;
230		hp->ns = n;
231	+	hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0;
232	+	memset(hp->sa, 0, sizeof(AMBSAMP)nn);
233	+	hp->sampOK = 0;
234		/* assign coefficient */
235	<	copycolor(hp->acoef, ac);
235	>	copycolor(hp->acoef, rcol);
236		d = 1.0/(n*n);
237		scalecolor(hp->acoef, d);
238	<	/* make tangent axes */
239	<	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
240	<	for (i = 0; i < 3; i++)
241	<	if (r->ron[i] < 0.6 && r->ron[i] > -0.6)
238	>	/* make tangent plane axes */
239	>	hp->uy[0] = 0.5 - frandom();
240	>	hp->uy[1] = 0.5 - frandom();
241	>	hp->uy[2] = 0.5 - frandom();
242	>	for (i = 3; i--; )
243	>	if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
244		break;
245	<	if (i >= 3)
246	<	error(CONSISTENCY, "bad ray direction in inithemi()");
245	>	if (i < 0)
246	>	error(CONSISTENCY, "bad ray direction in samp_hemi");
247		hp->uy[i] = 1.0;
248		VCROSS(hp->ux, hp->uy, r->ron);
249		normalize(hp->ux);
250		VCROSS(hp->uy, r->ron, hp->ux);
251	<	/* we're ready to sample */
252	<	return(hp);
251	>	/* sample divisions */
252	>	for (i = hp->ns; i--; )
253	>	for (j = hp->ns; j--; )
254	>	hp->sampOK += ambsample(hp, i, j, 0);
255	>	copycolor(rcol, hp->acol);
256	>	if (!hp->sampOK) { /* utter failure? */
257	>	free(hp);
258	>	return(NULL);
259	>	}
260	>	if (hp->sampOK < hp->ns*hp->ns) {
261	>	hp->sampOK = -1; / soft failure */
262	>	return(hp);
263	>	}
264	>	n = ambssamp*wt + 0.5;
265	>	if (n > 8) { /* perform super-sampling? */
266	>	ambsupersamp(hp, n);
267	>	copycolor(rcol, hp->acol);
268	>	}
269	>	return(hp); /* all is well */
270		}
271
272
273	<	static int
274	<	ambsample( /* sample an ambient direction */
275	<	AMBHEMI *hp,
91	<	int i,
92	<	int j
93	<	)
273	>	/* Return brightness of farthest ambient sample */
274	>	static double
275	>	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
276		{
277	<	struct s_ambsamp *ap = &ambsamp(hp,i,j);
278	<	RAY ar;
279	<	int hlist[3];
280	<	double spt[2], zd;
99	<	int ii;
100	<	/* ambient coefficient for weight */
101	<	if (ambacc > FTINY)
102	<	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
103	<	else
104	<	copycolor(ar.rcoef, hp->acoef);
105	<	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) {
106	<	setcolor(ap->v, 0., 0., 0.);
107	<	VCOPY(ap->p, hp->rp->rop);
108	<	return(0); /* no sample taken */
277	>	if (hp->sa[n1].d <= hp->sa[n2].d) {
278	>	if (hp->sa[n1].d <= hp->sa[n3].d)
279	>	return(colval(hp->sa[n1].v,CIEY));
280	>	return(colval(hp->sa[n3].v,CIEY));
281		}
282	<	if (ambacc > FTINY) {
283	<	multcolor(ar.rcoef, hp->acoef);
284	<	scalecolor(ar.rcoef, 1./AVGREFL);
113	<	}
114	<	/* generate hemispherical sample */
115	<	SDsquare2disk(spt, (i+.1+.8*frandom())/hp->ns,
116	<	(j+.1+.8*frandom())/hp->ns);
117	<	zd = sqrt(1. - spt[0]spt[0] - spt[1]spt[1]);
118	<	for (ii = 3; ii--; )
119	<	ar.rdir[ii] = spt[0]*hp->ux[ii] +
120	<	spt[1]*hp->uy[ii] +
121	<	zd*hp->rp->ron[ii];
122	<	checknorm(ar.rdir);
123	<	dimlist[ndims++] = i*hp->ns + j + 90171;
124	<	rayvalue(&ar); /* evaluate ray */
125	<	ndims--;
126	<	multcolor(ar.rcol, ar.rcoef); /* apply coefficient */
127	<	copycolor(ap->v, ar.rcol);
128	<	if (ar.rt > 20.0maxarad) / limit vertex distance */
129	<	ar.rt = 20.0*maxarad;
130	<	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
131	<	return(1);
282	>	if (hp->sa[n2].d <= hp->sa[n3].d)
283	>	return(colval(hp->sa[n2].v,CIEY));
284	>	return(colval(hp->sa[n3].v,CIEY));
285		}
286
287
288		/* Compute vectors and coefficients for Hessian/gradient calcs */
289		static void
290	<	comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop)
290	>	comp_fftri(FFTRI ftp, AMBHEMI hp, const int n0, const int n1)
291		{
292	<	FVECT v1;
293	<	double dot_e, dot_er, dot_r, dot_r1;
292	>	double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
293	>	int ii;
294
295	<	VSUB(ftp->r_i, ap0, rop);
296	<	VSUB(ftp->r_i1, ap1, rop);
297	<	VSUB(ftp->e_i, ap1, ap0);
298	<	VCROSS(v1, ftp->e_i, ftp->r_i);
299	<	ftp->nf = 1.0/DOT(v1,v1);
147	<	VCROSS(v1, ftp->r_i, ftp->r_i1);
148	<	ftp->I1 = sqrt(DOT(v1,v1)*ftp->nf);
295	>	VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop);
296	>	VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop);
297	>	VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p);
298	>	VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
299	>	rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
300		dot_e = DOT(ftp->e_i,ftp->e_i);
301		dot_er = DOT(ftp->e_i, ftp->r_i);
302	<	dot_r = DOT(ftp->r_i,ftp->r_i);
303	<	dot_r1 = DOT(ftp->r_i1,ftp->r_i1);
304	<	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r +
305	<	dot_eftp->I1 )0.5*ftp->nf;
306	<	ftp->J2 = 0.25ftp->nf( 1.0/dot_r - 1.0/dot_r1 ) -
307	<	dot_er/dot_e*ftp->I2;
302	>	rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i);
303	>	rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
304	>	ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_rrdot_r1) )
305	>	sqrt( rdot_cp );
306	>	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)rdot_r1 - dot_errdot_r +
307	>	dot_eftp->I1 )0.5*rdot_cp;
308	>	J2 = ( 0.5(rdot_r - rdot_r1) - dot_erftp->I2 ) / dot_e;
309	>	for (ii = 3; ii--; )
310	>	ftp->rI2_eJ2[ii] = ftp->I2ftp->r_i[ii] + J2ftp->e_i[ii];
311		}
312
313
314	<	/* Compose matrix from two vectors */
314	>	/* Compose 3x3 matrix from two vectors */
315		static void
316		compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
317		{
#	Line 174 \| Line 328 \| compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
328		static void
329		comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
330		{
331	<	FVECT v1, v2;
331	>	FVECT ncp;
332		FVECT m1[3], m2[3], m3[3], m4[3];
333		double d1, d2, d3, d4;
334		double I3, J3, K3;
#	Line 184 \| Line 338 \| *comp_hessian(FVECT hess[3], FFTRI ftp, FVECT nrm)**
338		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
339		d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
340		d4 = DOT(ftp->e_i, ftp->r_i);
341	<	I3 = 0.25ftp->nf( DOT(ftp->e_i, ftp->r_i1)d2d2 - d4d1d1 +
342	<	3.0ftp->I2d3 );
341	>	I3 = ( DOT(ftp->e_i, ftp->r_i1)d2d2 - d4d1d1 + 3.0/d3*ftp->I2 )
342	>	/ ( 4.0*DOT(ftp->rcp,ftp->rcp) );
343		J3 = 0.25d3(d1d1 - d2d2) - d4d3I3;
344		K3 = d3(ftp->I2 - I3/d1 - 2.0d4*J3);
345		/* intermediate matrices */
346	<	VCROSS(v1, nrm, ftp->e_i);
347	<	for (j = 3; j--; )
194	<	v2[i] = ftp->I2ftp->r_i[j] + ftp->J2ftp->e_i[j];
195	<	compose_matrix(m1, v1, v2);
346	>	VCROSS(ncp, nrm, ftp->e_i);
347	>	compose_matrix(m1, ncp, ftp->rI2_eJ2);
348		compose_matrix(m2, ftp->r_i, ftp->r_i);
349		compose_matrix(m3, ftp->e_i, ftp->e_i);
350		compose_matrix(m4, ftp->r_i, ftp->e_i);
351	<	VCROSS(v1, ftp->r_i, ftp->e_i);
200	<	d1 = DOT(nrm, v1);
351	>	d1 = DOT(nrm, ftp->rcp);
352		d2 = -d1*ftp->I2;
353		d1 *= 2.0;
354		for (i = 3; i--; ) /* final matrix sum */
#	Line 227 \| Line 378 \| rev_hessian(FVECT hess[3])
378		/* Add to radiometric Hessian from the given triangle */
379		static void
380		add2hessian(FVECT hess[3], FVECT ehess1[3],
381	<	FVECT ehess2[3], FVECT ehess3[3], COLORV v)
381	>	FVECT ehess2[3], FVECT ehess3[3], double v)
382		{
383		int i, j;
384
#	Line 241 \| Line 392 \| add2hessian(FVECT hess[3], FVECT ehess1[3],
392		static void
393		comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
394		{
395	<	FVECT vcp;
395	>	FVECT ncp;
396		double f1;
397		int i;
398
399	<	VCROSS(vcp, ftp->r_i, ftp->r_i1);
400	<	f1 = 2.0*DOT(nrm, vcp);
250	<	VCROSS(vcp, nrm, ftp->e_i);
399	>	f1 = 2.0*DOT(nrm, ftp->rcp);
400	>	VCROSS(ncp, nrm, ftp->e_i);
401		for (i = 3; i--; )
402	<	grad[i] = (0.5/PI)( ftp->I1vcp[i] +
253	<	f1(ftp->I2ftp->r_i[i] + ftp->J2*ftp->e_i[i]) );
402	>	grad[i] = (0.5/PI)( ftp->I1ncp[i] + f1*ftp->rI2_eJ2[i] );
403		}
404
405
#	Line 266 \| Line 415 \| rev_gradient(FVECT grad)
415
416		/* Add to displacement gradient from the given triangle */
417		static void
418	<	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
418	>	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
419		{
420		int i;
421
#	Line 275 \| Line 424 \| add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
424		}
425
426
278	–	/* Return brightness of furthest ambient sample */
279	–	static COLORV
280	–	back_ambval(struct s_ambsamp ap1, struct s_ambsamp ap2,
281	–	struct s_ambsamp *ap3, FVECT orig)
282	–	{
283	–	COLORV vback;
284	–	FVECT vec;
285	–	double d2, d2best;
286	–
287	–	VSUB(vec, ap1->p, orig);
288	–	d2best = DOT(vec,vec);
289	–	vback = ap1->v[CIEY];
290	–	VSUB(vec, ap2->p, orig);
291	–	d2 = DOT(vec,vec);
292	–	if (d2 > d2best) {
293	–	d2best = d2;
294	–	vback = ap2->v[CIEY];
295	–	}
296	–	VSUB(vec, ap3->p, orig);
297	–	d2 = DOT(vec,vec);
298	–	if (d2 > d2best)
299	–	return(ap3->v[CIEY]);
300	–	return(vback);
301	–	}
302	–
303	–
427		/* Compute anisotropic radii and eigenvector directions */
428	<	static int
428	>	static void
429		eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
430		{
431		double hess2[2][2];
#	Line 318 \| Line 441 \| eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
441		hess2[0][1] = DOT(uv[0], b);
442		hess2[1][0] = DOT(uv[1], a);
443		hess2[1][1] = DOT(uv[1], b);
444	<	/* compute eigenvalues */
445	<	if (quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
446	<	hess2[0][0]hess2[1][1]-hess2[0][1]hess2[1][0]) != 2 \|\|
447	<	(evalue[0] = fabs(evalue[0])) <= FTINYFTINYFTINY \|\|
448	<	(evalue[1] = fabs(evalue[1])) <= FTINYFTINYFTINY)
449	<	error(INTERNAL, "bad eigenvalue calculation");
450	<
444	>	/* compute eigenvalue(s) */
445	>	i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
446	>	hess2[0][0]hess2[1][1]-hess2[0][1]hess2[1][0]);
447	>	if (i == 1) /* double-root (circle) */
448	>	evalue[1] = evalue[0];
449	>	if (!i \|\| ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) \|
450	>	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) {
451	>	ra[0] = ra[1] = maxarad;
452	>	return;
453	>	}
454		if (evalue[0] > evalue[1]) {
455	<	ra[0] = 1.0/sqrt(sqrt(evalue[0]));
456	<	ra[1] = 1.0/sqrt(sqrt(evalue[1]));
455	>	ra[0] = sqrt(sqrt(4.0/evalue[0]));
456	>	ra[1] = sqrt(sqrt(4.0/evalue[1]));
457		slope1 = evalue[1];
458		} else {
459	<	ra[0] = 1.0/sqrt(sqrt(evalue[1]));
460	<	ra[1] = 1.0/sqrt(sqrt(evalue[0]));
459	>	ra[0] = sqrt(sqrt(4.0/evalue[1]));
460	>	ra[1] = sqrt(sqrt(4.0/evalue[0]));
461		slope1 = evalue[0];
462		}
463		/* compute unit eigenvectors */
#	Line 352 \| Line 478 \| static void
478		ambHessian( /* anisotropic radii & pos. gradient */
479		AMBHEMI *hp,
480		FVECT uv[2], /* returned */
481	<	float ra[2], /* returned */
482	<	float pg[2] /* returned */
481	>	float ra[2], /* returned (optional) */
482	>	float pg[2] /* returned (optional) */
483		)
484		{
485		static char memerrmsg[] = "out of memory in ambHessian()";
#	Line 368 \| Line 494 \| *ambHessian( / anisotropic radii & pos. gradient /*
494		VCOPY(uv[1], hp->uy);
495		/* clock-wise vertex traversal from sample POV */
496		if (ra != NULL) { /* initialize Hessian row buffer */
497	<	hessrow = (FVECT ()[3])malloc(sizeof(FVECT)3*hp->ns);
497	>	hessrow = (FVECT ()[3])malloc(sizeof(FVECT)3*(hp->ns-1));
498		if (hessrow == NULL)
499		error(SYSTEM, memerrmsg);
500		memset(hessian, 0, sizeof(hessian));
501		} else if (pg == NULL) /* bogus call? */
502		return;
503		if (pg != NULL) { /* initialize form factor row buffer */
504	<	gradrow = (FVECT )malloc(sizeof(FVECT)hp->ns);
504	>	gradrow = (FVECT )malloc(sizeof(FVECT)(hp->ns-1));
505		if (gradrow == NULL)
506		error(SYSTEM, memerrmsg);
507		memset(gradient, 0, sizeof(gradient));
508		}
509		/* compute first row of edges */
510		for (j = 0; j < hp->ns-1; j++) {
511	<	comp_fftri(&fftr, ambsamp(hp,0,j).p,
386	<	ambsamp(hp,0,j+1).p, hp->rp->rop);
511	>	comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1));
512		if (hessrow != NULL)
513		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
514		if (gradrow != NULL)
#	Line 393 \| Line 518 \| *ambHessian( / anisotropic radii & pos. gradient /*
518		for (i = 0; i < hp->ns-1; i++) {
519		FVECT hesscol[3]; /* compute first vertical edge */
520		FVECT gradcol;
521	<	comp_fftri(&fftr, ambsamp(hp,i,0).p,
397	<	ambsamp(hp,i+1,0).p, hp->rp->rop);
521	>	comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0));
522		if (hessrow != NULL)
523		comp_hessian(hesscol, &fftr, hp->rp->ron);
524		if (gradrow != NULL)
#	Line 402 \| Line 526 \| *ambHessian( / anisotropic radii & pos. gradient /*
526		for (j = 0; j < hp->ns-1; j++) {
527		FVECT hessdia[3]; /* compute triangle contributions */
528		FVECT graddia;
529	<	COLORV backg;
530	<	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
531	<	&ambsamp(hp,i+1,j), hp->rp->rop);
529	>	double backg;
530	>	backg = back_ambval(hp, AI(hp,i,j),
531	>	AI(hp,i,j+1), AI(hp,i+1,j));
532		/* diagonal (inner) edge */
533	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
410	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
533	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j));
534		if (hessrow != NULL) {
535		comp_hessian(hessdia, &fftr, hp->rp->ron);
536		rev_hessian(hesscol);
537		add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
538		}
539	<	if (gradient != NULL) {
539	>	if (gradrow != NULL) {
540		comp_gradient(graddia, &fftr, hp->rp->ron);
541		rev_gradient(gradcol);
542		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
543		}
544		/* initialize edge in next row */
545	<	comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
423	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
545	>	comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j));
546		if (hessrow != NULL)
547		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
548		if (gradrow != NULL)
549		comp_gradient(gradrow[j], &fftr, hp->rp->ron);
550		/* new column edge & paired triangle */
551	<	backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
552	<	&ambsamp(hp,i+1,j), hp->rp->rop);
553	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
432	<	hp->rp->rop);
551	>	backg = back_ambval(hp, AI(hp,i+1,j+1),
552	>	AI(hp,i+1,j), AI(hp,i,j+1));
553	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1));
554		if (hessrow != NULL) {
555		comp_hessian(hesscol, &fftr, hp->rp->ron);
556		rev_hessian(hessdia);
#	Line 452 \| Line 573 \| *ambHessian( / anisotropic radii & pos. gradient /*
573
574		if (ra != NULL) /* extract eigenvectors & radii */
575		eigenvectors(uv, ra, hessian);
576	<	if (pg != NULL) { /* project position gradient */
576	>	if (pg != NULL) { /* tangential position gradient */
577		pg[0] = DOT(gradient, uv[0]);
578		pg[1] = DOT(gradient, uv[1]);
579		}
#	Line 463 \| Line 584 \| *ambHessian( / anisotropic radii & pos. gradient /*
584		static void
585		ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
586		{
587	<	struct s_ambsamp *ap;
588	<	int n;
587	>	AMBSAMP *ap;
588	>	double dgsum[2];
589	>	int n;
590	>	FVECT vd;
591	>	double gfact;
592
593	<	dg[0] = dg[1] = 0;
593	>	dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */
594		for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
471	–	FVECT vd;
472	–	double gfact;
595		/* use vector for azimuth + 90deg */
596		VSUB(vd, ap->p, hp->rp->rop);
597	<	/* brightness with tangent factor */
598	<	gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd);
597	>	/* brightness over cosine factor */
598	>	gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
599		/* sine = proj_radius/vd_length */
600	<	dg[0] -= DOT(uv[1], vd) * gfact ;
601	<	dg[1] += DOT(uv[0], vd) * gfact;
600	>	dgsum[0] -= DOT(uv[1], vd) * gfact;
601	>	dgsum[1] += DOT(uv[0], vd) * gfact;
602		}
603	+	dg[0] = dgsum[0] / (hp->ns*hp->ns);
604	+	dg[1] = dgsum[1] / (hp->ns*hp->ns);
605		}
606
607
608	+	/* Compute potential light leak direction flags for cache value */
609	+	static uint32
610	+	ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
611	+	{
612	+	const double max_d = 1.0/(minarad*ambacc + 0.001);
613	+	const double ang_res = 0.5*PI/(hp->ns-1);
614	+	const double ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY));
615	+	double avg_d = 0;
616	+	uint32 flgs = 0;
617	+	FVECT vec;
618	+	double u, v;
619	+	double ang, a1;
620	+	int i, j;
621	+	/* don't bother for a few samples */
622	+	if (hp->ns < 12)
623	+	return(0);
624	+	/* check distances overhead */
625	+	for (i = hp->ns*3/4; i-- > hp->ns>>2; )
626	+	for (j = hp->ns*3/4; j-- > hp->ns>>2; )
627	+	avg_d += ambsam(hp,i,j).d;
628	+	avg_d = 4.0/(hp->nshp->ns);
629	+	if (avg_dr0 >= 1.0) / ceiling too low for corral? */
630	+	return(0);
631	+	if (avg_d >= max_d) /* insurance */
632	+	return(0);
633	+	/* else circle around perimeter */
634	+	for (i = 0; i < hp->ns; i++)
635	+	for (j = 0; j < hp->ns; j += !i\|(i==hp->ns-1) ? 1 : hp->ns-1) {
636	+	AMBSAMP *ap = &ambsam(hp,i,j);
637	+	if ((ap->d <= FTINY) \| (ap->d >= max_d))
638	+	continue; /* too far or too near */
639	+	VSUB(vec, ap->p, hp->rp->rop);
640	+	u = DOT(vec, uv[0]);
641	+	v = DOT(vec, uv[1]);
642	+	if ((r0r0uu + r1r1vv) * ap->dap->d <= uu + v*v)
643	+	continue; /* occluder outside ellipse */
644	+	ang = atan2a(v, u); /* else set direction flags */
645	+	for (a1 = ang-.5ang_res; a1 <= ang+.5ang_res; a1 += ang_step)
646	+	flgs \|= 1L<<(int)(16/PI(a1 + 2.PI*(a1 < 0)));
647	+	}
648	+	/* add low-angle incident (< 20deg) */
649	+	if (fabs(hp->rp->rod) <= 0.342) {
650	+	u = -DOT(hp->rp->rdir, uv[0]);
651	+	v = -DOT(hp->rp->rdir, uv[1]);
652	+	if ((r0r0uu + r1r1vv) > hp->rp->rot*hp->rp->rot) {
653	+	ang = atan2a(v, u);
654	+	ang += 2.PI(ang < 0);
655	+	ang *= 16/PI;
656	+	if ((ang < .5) \| (ang >= 31.5))
657	+	flgs \|= 0x80000001;
658	+	else
659	+	flgs \|= 3L<<(int)(ang-.5);
660	+	}
661	+	}
662	+	return(flgs);
663	+	}
664	+
665	+
666		int
667		doambient( /* compute ambient component */
668		COLOR rcol, /* input/output color */
#	Line 489 \| Line 671 \| *doambient( / compute ambient component /*
671		FVECT uv[2], /* returned (optional) */
672		float ra[2], /* returned (optional) */
673		float pg[2], /* returned (optional) */
674	<	float dg[2] /* returned (optional) */
674	>	float dg[2], /* returned (optional) */
675	>	uint32 crlp / returned (optional) */
676		)
677		{
678	<	int cnt = 0;
679	<	FVECT my_uv[2];
680	<	AMBHEMI *hp;
681	<	double d, acol[3];
682	<	struct s_ambsamp *ap;
683	<	int i, j;
501	<	/* initialize */
502	<	if ((hp = inithemi(rcol, r, wt)) == NULL)
503	<	return(0);
678	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
679	>	FVECT my_uv[2];
680	>	double d, K;
681	>	AMBSAMP *ap;
682	>	int i;
683	>	/* clear return values */
684		if (uv != NULL)
685		memset(uv, 0, sizeof(FVECT)*2);
686		if (ra != NULL)
#	Line 509 \| Line 689 \| *doambient( / compute ambient component /*
689		pg[0] = pg[1] = 0.0;
690		if (dg != NULL)
691		dg[0] = dg[1] = 0.0;
692	<	/* sample the hemisphere */
693	<	acol[0] = acol[1] = acol[2] = 0.0;
694	<	for (i = hp->ns; i--; )
695	<	for (j = hp->ns; j--; )
696	<	if (ambsample(hp, i, j)) {
697	<	ap = &ambsamp(hp,i,j);
698	<	addcolor(acol, ap->v);
699	<	++cnt;
700	<	}
521	<	if (!cnt) {
522	<	setcolor(rcol, 0.0, 0.0, 0.0);
523	<	free(hp);
524	<	return(0); /* no valid samples */
692	>	if (crlp != NULL)
693	>	*crlp = 0;
694	>	if (hp == NULL) /* sampling falure? */
695	>	return(0);
696	>
697	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) \|\|
698	>	(hp->sampOK < 0) \| (hp->ns < 4)) {
699	>	free(hp); /* Hessian not requested/possible */
700	>	return(-1); /* value-only return value */
701		}
702	<	d = 1.0 / cnt; /* final indirect irradiance/PI */
703	<	acol[0] = d; acol[1] = d; acol[2] *= d;
704	<	copycolor(rcol, acol);
705	<	if (cnt < hp->nshp->ns \|\| / incomplete sampling? */
706	<	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
707	<	free(hp);
708	<	return(-1); /* no radius or gradient calc. */
702	>	if ((d = bright(rcol)) > FTINY) { /* normalize Y values */
703	>	d = 0.99(hp->nshp->ns)/d;
704	>	K = 0.01;
705	>	} else { /* or fall back on geometric Hessian */
706	>	K = 1.0;
707	>	pg = NULL;
708	>	dg = NULL;
709	>	crlp = NULL;
710		}
711	<	d = 0.01 * bright(rcol); /* add in 1% before Hessian comp. */
535	<	if (d < FTINY) d = FTINY;
536	<	ap = hp->sa; /* using Y channel from here on... */
711	>	ap = hp->sa; /* relative Y channel from here on... */
712		for (i = hp->ns*hp->ns; i--; ap++)
713	<	colval(ap->v,CIEY) = bright(ap->v) + d;
713	>	colval(ap->v,CIEY) = bright(ap->v)*d + K;
714
715		if (uv == NULL) /* make sure we have axis pointers */
716		uv = my_uv;
717		/* compute radii & pos. gradient */
718		ambHessian(hp, uv, ra, pg);
719	+
720		if (dg != NULL) /* compute direction gradient */
721		ambdirgrad(hp, uv, dg);
722	<	if (ra != NULL) { /* adjust/clamp radii */
723	<	d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt));
724	<	if ((ra[0] *= d) > maxarad)
725	<	ra[0] = maxarad;
722	>
723	>	if (ra != NULL) { /* scale/clamp radii */
724	>	if (pg != NULL) {
725	>	if (ra[0]*(d = fabs(pg[0])) > 1.0)
726	>	ra[0] = 1.0/d;
727	>	if (ra[1]*(d = fabs(pg[1])) > 1.0)
728	>	ra[1] = 1.0/d;
729	>	if (ra[0] > ra[1])
730	>	ra[0] = ra[1];
731	>	}
732	>	if (ra[0] < minarad) {
733	>	ra[0] = minarad;
734	>	if (ra[1] < minarad)
735	>	ra[1] = minarad;
736	>	}
737	>	ra[0] *= d = 1.0/sqrt(wt);
738		if ((ra[1] = d) > 2.0ra[0])
739		ra[1] = 2.0*ra[0];
740	+	if (ra[1] > maxarad) {
741	+	ra[1] = maxarad;
742	+	if (ra[0] > maxarad)
743	+	ra[0] = maxarad;
744	+	}
745	+	/* flag encroached directions */
746	+	if ((wt >= 0.89*AVGREFL) & (crlp != NULL))
747	+	crlp = ambcorral(hp, uv, ra[0]ambacc, ra[1]*ambacc);
748	+	if (pg != NULL) { /* cap gradient if necessary */
749	+	d = pg[0]pg[0]ra[0]ra[0] + pg[1]pg[1]ra[1]ra[1];
750	+	if (d > 1.0) {
751	+	d = 1.0/sqrt(d);
752	+	pg[0] *= d;
753	+	pg[1] *= d;
754	+	}
755	+	}
756		}
757		free(hp); /* clean up and return */
758		return(1);

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing ray/src/rt/ambcomp.c (file contents): Revision 2.27 by greg, Sat Apr 19 02:39:44 2014 UTC vs. Revision 2.64 by greg, Tue Aug 19 15:04:40 2014 UTC

Diff Legend

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.27 by greg, Sat Apr 19 02:39:44 2014 UTC vs.
Revision 2.64 by greg, Tue Aug 19 15:04:40 2014 UTC