[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.71 by greg, Tue Mar 15 20:35:41 2016 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	+	float *ep;
184	+	int i, j, n, nss;
185	+
186	+	if (earr == NULL) /* just skip calc. if no memory */
187	+	return;
188	+	/* accumulate estimated variances */
189	+	for (ep = earr + hp->ns*hp->ns; ep > earr; )
190	+	e2rem += *--ep;
191	+	ep = earr; /* perform super-sampling */
192	+	for (ap = hp->sa, i = 0; i < hp->ns; i++)
193	+	for (j = 0; j < hp->ns; j++, ap++) {
194	+	if (e2rem <= FTINY)
195	+	goto done; /* nothing left to do */
196	+	nss = ep/e2remcnt + frandom();
197	+	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
198	+	--cnt;
199	+	e2rem -= ep++; / update remainder */
200	+	}
201	+	done:
202	+	free(earr);
203	+	}
204	+
205	+
206		static AMBHEMI *
207	<	inithemi( /* initialize sampling hemisphere */
208	<	COLOR ac,
207	>	samp_hemi( /* sample indirect hemisphere */
208	>	COLOR rcol,
209		RAY *r,
210		double wt
211		)
212		{
213		AMBHEMI *hp;
214		double d;
215	<	int n, i;
215	>	int n, i, j;
216		/* set number of divisions */
217		if (ambacc <= FTINY &&
218	<	wt > (d = 0.8intens(ac)r->rweight/(ambdiv*minweight)))
218	>	wt > (d = 0.8intens(rcol)r->rweight/(ambdiv*minweight)))
219		wt = d; /* avoid ray termination */
220		n = sqrt(ambdiv * wt) + 0.5;
221		i = 1 + 5(ambacc > FTINY); / minimum number of samples */
222		if (n < i)
223		n = i;
224		/* allocate sampling array */
225	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
63	<	sizeof(struct s_ambsamp)(nn - 1));
225	>	hp = (AMBHEMI )malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)(n*n - 1));
226		if (hp == NULL)
227	<	return(NULL);
227	>	error(SYSTEM, "out of memory in samp_hemi");
228		hp->rp = r;
229		hp->ns = n;
230	+	hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0;
231	+	memset(hp->sa, 0, sizeof(AMBSAMP)nn);
232	+	hp->sampOK = 0;
233		/* assign coefficient */
234	<	copycolor(hp->acoef, ac);
234	>	copycolor(hp->acoef, rcol);
235		d = 1.0/(n*n);
236		scalecolor(hp->acoef, d);
237	<	/* make tangent axes */
238	<	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
239	<	for (i = 0; i < 3; i++)
75	<	if (r->ron[i] < 0.6 && r->ron[i] > -0.6)
76	<	break;
77	<	if (i >= 3)
78	<	error(CONSISTENCY, "bad ray direction in inithemi()");
79	<	hp->uy[i] = 1.0;
80	<	VCROSS(hp->ux, hp->uy, r->ron);
81	<	normalize(hp->ux);
237	>	/* make tangent plane axes */
238	>	if (!getperpendicular(hp->ux, r->ron, 1))
239	>	error(CONSISTENCY, "bad ray direction in samp_hemi");
240		VCROSS(hp->uy, r->ron, hp->ux);
241	<	/* we're ready to sample */
242	<	return(hp);
241	>	/* sample divisions */
242	>	for (i = hp->ns; i--; )
243	>	for (j = hp->ns; j--; )
244	>	hp->sampOK += ambsample(hp, i, j, 0);
245	>	copycolor(rcol, hp->acol);
246	>	if (!hp->sampOK) { /* utter failure? */
247	>	free(hp);
248	>	return(NULL);
249	>	}
250	>	if (hp->sampOK < hp->ns*hp->ns) {
251	>	hp->sampOK = -1; / soft failure */
252	>	return(hp);
253	>	}
254	>	n = ambssamp*wt + 0.5;
255	>	if (n > 8) { /* perform super-sampling? */
256	>	ambsupersamp(hp, n);
257	>	copycolor(rcol, hp->acol);
258	>	}
259	>	return(hp); /* all is well */
260		}
261
262
263	<	static int
264	<	ambsample( /* sample an ambient direction */
265	<	AMBHEMI *hp,
91	<	int i,
92	<	int j
93	<	)
263	>	/* Return brightness of farthest ambient sample */
264	>	static double
265	>	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
266		{
267	<	struct s_ambsamp *ap = &ambsamp(hp,i,j);
268	<	RAY ar;
269	<	int hlist[3];
270	<	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 */
267	>	if (hp->sa[n1].d <= hp->sa[n2].d) {
268	>	if (hp->sa[n1].d <= hp->sa[n3].d)
269	>	return(colval(hp->sa[n1].v,CIEY));
270	>	return(colval(hp->sa[n3].v,CIEY));
271		}
272	<	if (ambacc > FTINY) {
273	<	multcolor(ar.rcoef, hp->acoef);
274	<	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);
272	>	if (hp->sa[n2].d <= hp->sa[n3].d)
273	>	return(colval(hp->sa[n2].v,CIEY));
274	>	return(colval(hp->sa[n3].v,CIEY));
275		}
276
277
278		/* Compute vectors and coefficients for Hessian/gradient calcs */
279		static void
280	<	comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop)
280	>	comp_fftri(FFTRI ftp, AMBHEMI hp, const int n0, const int n1)
281		{
282	<	FVECT v1;
283	<	double dot_e, dot_er, dot_r, dot_r1;
282	>	double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
283	>	int ii;
284
285	<	VSUB(ftp->r_i, ap0, rop);
286	<	VSUB(ftp->r_i1, ap1, rop);
287	<	VSUB(ftp->e_i, ap1, ap0);
288	<	VCROSS(v1, ftp->e_i, ftp->r_i);
289	<	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);
285	>	VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop);
286	>	VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop);
287	>	VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p);
288	>	VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
289	>	rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
290		dot_e = DOT(ftp->e_i,ftp->e_i);
291		dot_er = DOT(ftp->e_i, ftp->r_i);
292	<	dot_r = DOT(ftp->r_i,ftp->r_i);
293	<	dot_r1 = DOT(ftp->r_i1,ftp->r_i1);
294	<	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r +
295	<	dot_eftp->I1 )0.5*ftp->nf;
296	<	ftp->J2 = 0.25ftp->nf( 1.0/dot_r - 1.0/dot_r1 ) -
297	<	dot_er/dot_e*ftp->I2;
292	>	rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i);
293	>	rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
294	>	ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_rrdot_r1) )
295	>	sqrt( rdot_cp );
296	>	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)rdot_r1 - dot_errdot_r +
297	>	dot_eftp->I1 )0.5*rdot_cp;
298	>	J2 = ( 0.5(rdot_r - rdot_r1) - dot_erftp->I2 ) / dot_e;
299	>	for (ii = 3; ii--; )
300	>	ftp->rI2_eJ2[ii] = ftp->I2ftp->r_i[ii] + J2ftp->e_i[ii];
301		}
302
303
304	<	/* Compose matrix from two vectors */
304	>	/* Compose 3x3 matrix from two vectors */
305		static void
306		compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
307		{
#	Line 174 \| Line 318 \| compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
318		static void
319		comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
320		{
321	<	FVECT v1, v2;
321	>	FVECT ncp;
322		FVECT m1[3], m2[3], m3[3], m4[3];
323		double d1, d2, d3, d4;
324		double I3, J3, K3;
#	Line 184 \| Line 328 \| *comp_hessian(FVECT hess[3], FFTRI ftp, FVECT nrm)**
328		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
329		d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
330		d4 = DOT(ftp->e_i, ftp->r_i);
331	<	I3 = 0.25ftp->nf( DOT(ftp->e_i, ftp->r_i1)d2d2 - d4d1d1 +
332	<	3.0ftp->I2d3 );
331	>	I3 = ( DOT(ftp->e_i, ftp->r_i1)d2d2 - d4d1d1 + 3.0/d3*ftp->I2 )
332	>	/ ( 4.0*DOT(ftp->rcp,ftp->rcp) );
333		J3 = 0.25d3(d1d1 - d2d2) - d4d3I3;
334		K3 = d3(ftp->I2 - I3/d1 - 2.0d4*J3);
335		/* intermediate matrices */
336	<	VCROSS(v1, nrm, ftp->e_i);
337	<	for (j = 3; j--; )
194	<	v2[i] = ftp->I2ftp->r_i[j] + ftp->J2ftp->e_i[j];
195	<	compose_matrix(m1, v1, v2);
336	>	VCROSS(ncp, nrm, ftp->e_i);
337	>	compose_matrix(m1, ncp, ftp->rI2_eJ2);
338		compose_matrix(m2, ftp->r_i, ftp->r_i);
339		compose_matrix(m3, ftp->e_i, ftp->e_i);
340		compose_matrix(m4, ftp->r_i, ftp->e_i);
341	<	VCROSS(v1, ftp->r_i, ftp->e_i);
200	<	d1 = DOT(nrm, v1);
341	>	d1 = DOT(nrm, ftp->rcp);
342		d2 = -d1*ftp->I2;
343		d1 *= 2.0;
344		for (i = 3; i--; ) /* final matrix sum */
#	Line 227 \| Line 368 \| rev_hessian(FVECT hess[3])
368		/* Add to radiometric Hessian from the given triangle */
369		static void
370		add2hessian(FVECT hess[3], FVECT ehess1[3],
371	<	FVECT ehess2[3], FVECT ehess3[3], COLORV v)
371	>	FVECT ehess2[3], FVECT ehess3[3], double v)
372		{
373		int i, j;
374
#	Line 241 \| Line 382 \| add2hessian(FVECT hess[3], FVECT ehess1[3],
382		static void
383		comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
384		{
385	<	FVECT vcp;
385	>	FVECT ncp;
386		double f1;
387		int i;
388
389	<	VCROSS(vcp, ftp->r_i, ftp->r_i1);
390	<	f1 = 2.0*DOT(nrm, vcp);
250	<	VCROSS(vcp, nrm, ftp->e_i);
389	>	f1 = 2.0*DOT(nrm, ftp->rcp);
390	>	VCROSS(ncp, nrm, ftp->e_i);
391		for (i = 3; i--; )
392	<	grad[i] = (0.5/PI)( ftp->I1vcp[i] +
253	<	f1(ftp->I2ftp->r_i[i] + ftp->J2*ftp->e_i[i]) );
392	>	grad[i] = (0.5/PI)( ftp->I1ncp[i] + f1*ftp->rI2_eJ2[i] );
393		}
394
395
#	Line 266 \| Line 405 \| rev_gradient(FVECT grad)
405
406		/* Add to displacement gradient from the given triangle */
407		static void
408	<	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
408	>	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
409		{
410		int i;
411
#	Line 275 \| Line 414 \| add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
414		}
415
416
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	–
417		/* Compute anisotropic radii and eigenvector directions */
418	<	static int
418	>	static void
419		eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
420		{
421		double hess2[2][2];
#	Line 318 \| Line 431 \| eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
431		hess2[0][1] = DOT(uv[0], b);
432		hess2[1][0] = DOT(uv[1], a);
433		hess2[1][1] = DOT(uv[1], b);
434	<	/* compute eigenvalues */
435	<	if (quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
436	<	hess2[0][0]hess2[1][1]-hess2[0][1]hess2[1][0]) != 2 \|\|
437	<	(evalue[0] = fabs(evalue[0])) <= FTINYFTINYFTINY \|\|
438	<	(evalue[1] = fabs(evalue[1])) <= FTINYFTINYFTINY)
439	<	error(INTERNAL, "bad eigenvalue calculation");
440	<
434	>	/* compute eigenvalue(s) */
435	>	i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
436	>	hess2[0][0]hess2[1][1]-hess2[0][1]hess2[1][0]);
437	>	if (i == 1) /* double-root (circle) */
438	>	evalue[1] = evalue[0];
439	>	if (!i \|\| ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) \|
440	>	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) {
441	>	ra[0] = ra[1] = maxarad;
442	>	return;
443	>	}
444		if (evalue[0] > evalue[1]) {
445	<	ra[0] = 1.0/sqrt(sqrt(evalue[0]));
446	<	ra[1] = 1.0/sqrt(sqrt(evalue[1]));
445	>	ra[0] = sqrt(sqrt(4.0/evalue[0]));
446	>	ra[1] = sqrt(sqrt(4.0/evalue[1]));
447		slope1 = evalue[1];
448		} else {
449	<	ra[0] = 1.0/sqrt(sqrt(evalue[1]));
450	<	ra[1] = 1.0/sqrt(sqrt(evalue[0]));
449	>	ra[0] = sqrt(sqrt(4.0/evalue[1]));
450	>	ra[1] = sqrt(sqrt(4.0/evalue[0]));
451		slope1 = evalue[0];
452		}
453		/* compute unit eigenvectors */
#	Line 352 \| Line 468 \| static void
468		ambHessian( /* anisotropic radii & pos. gradient */
469		AMBHEMI *hp,
470		FVECT uv[2], /* returned */
471	<	float ra[2], /* returned */
472	<	float pg[2] /* returned */
471	>	float ra[2], /* returned (optional) */
472	>	float pg[2] /* returned (optional) */
473		)
474		{
475		static char memerrmsg[] = "out of memory in ambHessian()";
#	Line 368 \| Line 484 \| *ambHessian( / anisotropic radii & pos. gradient /*
484		VCOPY(uv[1], hp->uy);
485		/* clock-wise vertex traversal from sample POV */
486		if (ra != NULL) { /* initialize Hessian row buffer */
487	<	hessrow = (FVECT ()[3])malloc(sizeof(FVECT)3*hp->ns);
487	>	hessrow = (FVECT ()[3])malloc(sizeof(FVECT)3*(hp->ns-1));
488		if (hessrow == NULL)
489		error(SYSTEM, memerrmsg);
490		memset(hessian, 0, sizeof(hessian));
491		} else if (pg == NULL) /* bogus call? */
492		return;
493		if (pg != NULL) { /* initialize form factor row buffer */
494	<	gradrow = (FVECT )malloc(sizeof(FVECT)hp->ns);
494	>	gradrow = (FVECT )malloc(sizeof(FVECT)(hp->ns-1));
495		if (gradrow == NULL)
496		error(SYSTEM, memerrmsg);
497		memset(gradient, 0, sizeof(gradient));
498		}
499		/* compute first row of edges */
500		for (j = 0; j < hp->ns-1; j++) {
501	<	comp_fftri(&fftr, ambsamp(hp,0,j).p,
386	<	ambsamp(hp,0,j+1).p, hp->rp->rop);
501	>	comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1));
502		if (hessrow != NULL)
503		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
504		if (gradrow != NULL)
#	Line 393 \| Line 508 \| *ambHessian( / anisotropic radii & pos. gradient /*
508		for (i = 0; i < hp->ns-1; i++) {
509		FVECT hesscol[3]; /* compute first vertical edge */
510		FVECT gradcol;
511	<	comp_fftri(&fftr, ambsamp(hp,i,0).p,
397	<	ambsamp(hp,i+1,0).p, hp->rp->rop);
511	>	comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0));
512		if (hessrow != NULL)
513		comp_hessian(hesscol, &fftr, hp->rp->ron);
514		if (gradrow != NULL)
#	Line 402 \| Line 516 \| *ambHessian( / anisotropic radii & pos. gradient /*
516		for (j = 0; j < hp->ns-1; j++) {
517		FVECT hessdia[3]; /* compute triangle contributions */
518		FVECT graddia;
519	<	COLORV backg;
520	<	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
521	<	&ambsamp(hp,i+1,j), hp->rp->rop);
519	>	double backg;
520	>	backg = back_ambval(hp, AI(hp,i,j),
521	>	AI(hp,i,j+1), AI(hp,i+1,j));
522		/* diagonal (inner) edge */
523	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
410	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
523	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j));
524		if (hessrow != NULL) {
525		comp_hessian(hessdia, &fftr, hp->rp->ron);
526		rev_hessian(hesscol);
527		add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
528		}
529	<	if (gradient != NULL) {
529	>	if (gradrow != NULL) {
530		comp_gradient(graddia, &fftr, hp->rp->ron);
531		rev_gradient(gradcol);
532		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
533		}
534		/* initialize edge in next row */
535	<	comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
423	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
535	>	comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j));
536		if (hessrow != NULL)
537		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
538		if (gradrow != NULL)
539		comp_gradient(gradrow[j], &fftr, hp->rp->ron);
540		/* new column edge & paired triangle */
541	<	backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
542	<	&ambsamp(hp,i+1,j), hp->rp->rop);
543	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
432	<	hp->rp->rop);
541	>	backg = back_ambval(hp, AI(hp,i+1,j+1),
542	>	AI(hp,i+1,j), AI(hp,i,j+1));
543	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1));
544		if (hessrow != NULL) {
545		comp_hessian(hesscol, &fftr, hp->rp->ron);
546		rev_hessian(hessdia);
#	Line 452 \| Line 563 \| *ambHessian( / anisotropic radii & pos. gradient /*
563
564		if (ra != NULL) /* extract eigenvectors & radii */
565		eigenvectors(uv, ra, hessian);
566	<	if (pg != NULL) { /* project position gradient */
566	>	if (pg != NULL) { /* tangential position gradient */
567		pg[0] = DOT(gradient, uv[0]);
568		pg[1] = DOT(gradient, uv[1]);
569		}
#	Line 463 \| Line 574 \| *ambHessian( / anisotropic radii & pos. gradient /*
574		static void
575		ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
576		{
577	<	struct s_ambsamp *ap;
578	<	int n;
577	>	AMBSAMP *ap;
578	>	double dgsum[2];
579	>	int n;
580	>	FVECT vd;
581	>	double gfact;
582
583	<	dg[0] = dg[1] = 0;
583	>	dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */
584		for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
471	–	FVECT vd;
472	–	double gfact;
585		/* use vector for azimuth + 90deg */
586		VSUB(vd, ap->p, hp->rp->rop);
587	<	/* brightness with tangent factor */
588	<	gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd);
587	>	/* brightness over cosine factor */
588	>	gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
589		/* sine = proj_radius/vd_length */
590	<	dg[0] -= DOT(uv[1], vd) * gfact ;
591	<	dg[1] += DOT(uv[0], vd) * gfact;
590	>	dgsum[0] -= DOT(uv[1], vd) * gfact;
591	>	dgsum[1] += DOT(uv[0], vd) * gfact;
592		}
593	+	dg[0] = dgsum[0] / (hp->ns*hp->ns);
594	+	dg[1] = dgsum[1] / (hp->ns*hp->ns);
595		}
596
597
598	+	/* Compute potential light leak direction flags for cache value */
599	+	static uint32
600	+	ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
601	+	{
602	+	const double max_d = 1.0/(minarad*ambacc + 0.001);
603	+	const double ang_res = 0.5*PI/hp->ns;
604	+	const double ang_step = ang_res/((int)(16/PI*ang_res) + 1.01);
605	+	double avg_d = 0;
606	+	uint32 flgs = 0;
607	+	FVECT vec;
608	+	double u, v;
609	+	double ang, a1;
610	+	int i, j;
611	+	/* don't bother for a few samples */
612	+	if (hp->ns < 12)
613	+	return(0);
614	+	/* check distances overhead */
615	+	for (i = hp->ns*3/4; i-- > hp->ns>>2; )
616	+	for (j = hp->ns*3/4; j-- > hp->ns>>2; )
617	+	avg_d += ambsam(hp,i,j).d;
618	+	avg_d = 4.0/(hp->nshp->ns);
619	+	if (avg_dr0 >= 1.0) / ceiling too low for corral? */
620	+	return(0);
621	+	if (avg_d >= max_d) /* insurance */
622	+	return(0);
623	+	/* else circle around perimeter */
624	+	for (i = 0; i < hp->ns; i++)
625	+	for (j = 0; j < hp->ns; j += !i\|(i==hp->ns-1) ? 1 : hp->ns-1) {
626	+	AMBSAMP *ap = &ambsam(hp,i,j);
627	+	if ((ap->d <= FTINY) \| (ap->d >= max_d))
628	+	continue; /* too far or too near */
629	+	VSUB(vec, ap->p, hp->rp->rop);
630	+	u = DOT(vec, uv[0]);
631	+	v = DOT(vec, uv[1]);
632	+	if ((r0r0uu + r1r1vv) * ap->dap->d <= uu + v*v)
633	+	continue; /* occluder outside ellipse */
634	+	ang = atan2a(v, u); /* else set direction flags */
635	+	for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step)
636	+	flgs \|= 1L<<(int)(16/PI(a1 + 2.PI*(a1 < 0)));
637	+	}
638	+	/* add low-angle incident (< 20deg) */
639	+	if (fabs(hp->rp->rod) <= 0.342) {
640	+	u = -DOT(hp->rp->rdir, uv[0]);
641	+	v = -DOT(hp->rp->rdir, uv[1]);
642	+	if ((r0r0uu + r1r1vv) > hp->rp->rot*hp->rp->rot) {
643	+	ang = atan2a(v, u);
644	+	ang += 2.PI(ang < 0);
645	+	ang *= 16/PI;
646	+	if ((ang < .5) \| (ang >= 31.5))
647	+	flgs \|= 0x80000001;
648	+	else
649	+	flgs \|= 3L<<(int)(ang-.5);
650	+	}
651	+	}
652	+	return(flgs);
653	+	}
654	+
655	+
656		int
657		doambient( /* compute ambient component */
658		COLOR rcol, /* input/output color */
#	Line 489 \| Line 661 \| *doambient( / compute ambient component /*
661		FVECT uv[2], /* returned (optional) */
662		float ra[2], /* returned (optional) */
663		float pg[2], /* returned (optional) */
664	<	float dg[2] /* returned (optional) */
664	>	float dg[2], /* returned (optional) */
665	>	uint32 crlp / returned (optional) */
666		)
667		{
668	<	int cnt = 0;
669	<	FVECT my_uv[2];
670	<	AMBHEMI *hp;
671	<	double d, acol[3];
672	<	struct s_ambsamp *ap;
673	<	int i, j;
501	<	/* initialize */
502	<	if ((hp = inithemi(rcol, r, wt)) == NULL)
503	<	return(0);
668	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
669	>	FVECT my_uv[2];
670	>	double d, K;
671	>	AMBSAMP *ap;
672	>	int i;
673	>	/* clear return values */
674		if (uv != NULL)
675		memset(uv, 0, sizeof(FVECT)*2);
676		if (ra != NULL)
#	Line 509 \| Line 679 \| *doambient( / compute ambient component /*
679		pg[0] = pg[1] = 0.0;
680		if (dg != NULL)
681		dg[0] = dg[1] = 0.0;
682	<	/* sample the hemisphere */
683	<	acol[0] = acol[1] = acol[2] = 0.0;
684	<	for (i = hp->ns; i--; )
685	<	for (j = hp->ns; j--; )
686	<	if (ambsample(hp, i, j)) {
687	<	ap = &ambsamp(hp,i,j);
688	<	addcolor(acol, ap->v);
689	<	++cnt;
690	<	}
521	<	if (!cnt) {
522	<	setcolor(rcol, 0.0, 0.0, 0.0);
523	<	free(hp);
524	<	return(0); /* no valid samples */
682	>	if (crlp != NULL)
683	>	*crlp = 0;
684	>	if (hp == NULL) /* sampling falure? */
685	>	return(0);
686	>
687	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) \|\|
688	>	(hp->sampOK < 0) \| (hp->ns < 6)) {
689	>	free(hp); /* Hessian not requested/possible */
690	>	return(-1); /* value-only return value */
691		}
692	<	d = 1.0 / cnt; /* final indirect irradiance/PI */
693	<	acol[0] = d; acol[1] = d; acol[2] *= d;
694	<	copycolor(rcol, acol);
695	<	if (cnt < hp->nshp->ns \|\| / incomplete sampling? */
696	<	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
697	<	free(hp);
698	<	return(-1); /* no radius or gradient calc. */
692	>	if ((d = bright(rcol)) > FTINY) { /* normalize Y values */
693	>	d = 0.99(hp->nshp->ns)/d;
694	>	K = 0.01;
695	>	} else { /* or fall back on geometric Hessian */
696	>	K = 1.0;
697	>	pg = NULL;
698	>	dg = NULL;
699	>	crlp = NULL;
700		}
701	<	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... */
701	>	ap = hp->sa; /* relative Y channel from here on... */
702		for (i = hp->ns*hp->ns; i--; ap++)
703	<	colval(ap->v,CIEY) = bright(ap->v) + d;
703	>	colval(ap->v,CIEY) = bright(ap->v)*d + K;
704
705		if (uv == NULL) /* make sure we have axis pointers */
706		uv = my_uv;
707		/* compute radii & pos. gradient */
708		ambHessian(hp, uv, ra, pg);
709	+
710		if (dg != NULL) /* compute direction gradient */
711		ambdirgrad(hp, uv, dg);
712	<	if (ra != NULL) { /* adjust/clamp radii */
713	<	d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt));
714	<	if ((ra[0] *= d) > maxarad)
715	<	ra[0] = maxarad;
712	>
713	>	if (ra != NULL) { /* scale/clamp radii */
714	>	if (pg != NULL) {
715	>	if (ra[0]*(d = fabs(pg[0])) > 1.0)
716	>	ra[0] = 1.0/d;
717	>	if (ra[1]*(d = fabs(pg[1])) > 1.0)
718	>	ra[1] = 1.0/d;
719	>	if (ra[0] > ra[1])
720	>	ra[0] = ra[1];
721	>	}
722	>	if (ra[0] < minarad) {
723	>	ra[0] = minarad;
724	>	if (ra[1] < minarad)
725	>	ra[1] = minarad;
726	>	}
727	>	ra[0] *= d = 1.0/sqrt(wt);
728		if ((ra[1] = d) > 2.0ra[0])
729		ra[1] = 2.0*ra[0];
730	+	if (ra[1] > maxarad) {
731	+	ra[1] = maxarad;
732	+	if (ra[0] > maxarad)
733	+	ra[0] = maxarad;
734	+	}
735	+	/* flag encroached directions */
736	+	if (crlp != NULL)
737	+	crlp = ambcorral(hp, uv, ra[0]ambacc, ra[1]*ambacc);
738	+	if (pg != NULL) { /* cap gradient if necessary */
739	+	d = pg[0]pg[0]ra[0]ra[0] + pg[1]pg[1]ra[1]ra[1];
740	+	if (d > 1.0) {
741	+	d = 1.0/sqrt(d);
742	+	pg[0] *= d;
743	+	pg[1] *= d;
744	+	}
745	+	}
746		}
747		free(hp); /* clean up and return */
748		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.71 by greg, Tue Mar 15 20:35:41 2016 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.71 by greg, Tue Mar 15 20:35:41 2016 UTC