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.32 by greg, Thu Apr 24 17:36:43 2014 UTC vs.
Revision 2.76 by greg, Thu Jan 26 16:46:58 2017 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	<	FVECT   p;              /* 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, rI2_eJ2;
49	<	double  nf, I1, I2;
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 AMBHEMI *
54	<	inithemi(                       /* initialize sampling hemisphere */
55	<	COLOR   ac,
56	<	RAY     *r,
57	<	double  wt
53	>	static int
54	>	ambcollision(                           /* proposed direciton collides? */
55	>	AMBHEMI *hp,
56	>	int     i,
57	>	int     j,
58	>	FVECT   dv
59		)
60		{
61	<	AMBHEMI *hp;
62	<	double  d;
63	<	int     n, i;
64	<	/* set number of divisions */
65	<	if (ambacc <= FTINY &&
66	<	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
67	<	wt = d;                 /* avoid ray termination */
68	<	n = sqrt(ambdiv * wt) + 0.5;
69	<	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
70	<	if (n < i)
71	<	n = i;
72	<	/* allocate sampling array */
73	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
74	<	sizeof(struct s_ambsamp)*(n*n - 1));
75	<	if (hp == NULL)
76	<	return(NULL);
77	<	hp->rp = r;
78	<	hp->ns = n;
79	<	/* assign coefficient */
80	<	copycolor(hp->acoef, ac);
81	<	d = 1.0/(n*n);
82	<	scalecolor(hp->acoef, d);
83	<	/* make tangent plane axes */
84	<	hp->uy[0] = 0.1 - 0.2*frandom();
85	<	hp->uy[1] = 0.1 - 0.2*frandom();
86	<	hp->uy[2] = 0.1 - 0.2*frandom();
76	<	for (i = 0; i < 3; i++)
77	<	if (r->ron[i] < 0.6 && r->ron[i] > -0.6)
78	<	break;
79	<	if (i >= 3)
80	<	error(CONSISTENCY, "bad ray direction in inithemi()");
81	<	hp->uy[i] = 1.0;
82	<	VCROSS(hp->ux, hp->uy, r->ron);
83	<	normalize(hp->ux);
84	<	VCROSS(hp->uy, r->ron, hp->ux);
85	<	/* we're ready to sample */
86	<	return(hp);
61	>	double  cos_thresh;
62	>	int     ii, jj;
63	>	/* min. spacing = 1/4th division */
64	>	cos_thresh = (PI/4.)/(double)hp->ns;
65	>	cos_thresh = 1. - .5*cos_thresh*cos_thresh;
66	>	/* check existing neighbors */
67	>	for (ii = i-1; ii <= i+1; ii++) {
68	>	if (ii < 0) continue;
69	>	if (ii >= hp->ns) break;
70	>	for (jj = j-1; jj <= j+1; jj++) {
71	>	AMBSAMP *ap;
72	>	FVECT   avec;
73	>	double  dprod;
74	>	if (jj < 0) continue;
75	>	if (jj >= hp->ns) break;
76	>	if ((ii==i) & (jj==j)) continue;
77	>	ap = &ambsam(hp,ii,jj);
78	>	if (ap->d <= .5/FHUGE)
79	>	continue;       /* no one home */
80	>	VSUB(avec, ap->p, hp->rp->rop);
81	>	dprod = DOT(avec, dv);
82	>	if (dprod >= cos_thresh*VLEN(avec))
83	>	return(1);      /* collision */
84	>	}
85	>	}
86	>	return(0);                      /* nothing to worry about */
87		}
88
89
90	<	static struct s_ambsamp *
91	<	ambsample(                              /* sample an ambient direction */
90	>	static int
91	>	ambsample(                              /* initial ambient division sample */
92		AMBHEMI *hp,
93		int     i,
94	<	int     j
94	>	int     j,
95	>	int     n
96		)
97		{
98	<	struct s_ambsamp        *ap = &ambsamp(hp,i,j);
99	<	RAY                     ar;
100	<	double                  spt[2], zd;
101	<	int                     ii;
98	>	AMBSAMP *ap = &ambsam(hp,i,j);
99	>	RAY     ar;
100	>	int     hlist[3], ii;
101	>	double  spt[2], zd;
102	>	/* generate hemispherical sample */
103		/* ambient coefficient for weight */
104		if (ambacc > FTINY)
105		setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
106		else
107		copycolor(ar.rcoef, hp->acoef);
108		if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
109	<	goto badsample;
109	>	return(0);
110		if (ambacc > FTINY) {
111		multcolor(ar.rcoef, hp->acoef);
112		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 );
114	>	hlist[0] = hp->rp->rno;
115	>	hlist[1] = j;
116	>	hlist[2] = i;
117	>	multisamp(spt, 2, urand(ilhash(hlist,3)+n));
118	>	resample:
119	>	SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
120		zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
121		for (ii = 3; ii--; )
122		ar.rdir[ii] =   spt[0]*hp->ux[ii] +
123		spt[1]*hp->uy[ii] +
124		zd*hp->rp->ron[ii];
125		checknorm(ar.rdir);
126	<	dimlist[ndims++] = i*hp->ns + j + 90171;
126	>	/* avoid coincident samples */
127	>	if (!n && ambcollision(hp, i, j, ar.rdir)) {
128	>	spt[0] = frandom(); spt[1] = frandom();
129	>	goto resample;          /* reject this sample */
130	>	}
131	>	dimlist[ndims++] = AI(hp,i,j) + 90171;
132		rayvalue(&ar);                  /* evaluate ray */
133		ndims--;
134	<	if (ar.rt > 20.0*maxarad)       /* limit vertex distance */
135	<	ar.rt = 20.0*maxarad;
126	<	else if (ar.rt <= FTINY)        /* should never happen! */
127	<	goto badsample;
128	<	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
134	>	if (ar.rt <= FTINY)
135	>	return(0);              /* should never happen */
136		multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
137	<	copycolor(ap->v, ar.rcol);
138	<	return(ap);
139	<	badsample:
140	<	setcolor(ap->v, 0., 0., 0.);
141	<	VCOPY(ap->p, hp->rp->rop);
142	<	return(NULL);
137	>	if (ar.rt*ap->d < 1.0)          /* new/closer distance? */
138	>	ap->d = 1.0/ar.rt;
139	>	if (!n) {                       /* record first vertex & value */
140	>	if (ar.rt > 10.0*thescene.cusize + 1000.)
141	>	ar.rt = 10.0*thescene.cusize + 1000.;
142	>	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
143	>	copycolor(ap->v, ar.rcol);
144	>	} else {                        /* else update recorded value */
145	>	hp->acol[RED] -= colval(ap->v,RED);
146	>	hp->acol[GRN] -= colval(ap->v,GRN);
147	>	hp->acol[BLU] -= colval(ap->v,BLU);
148	>	zd = 1.0/(double)(n+1);
149	>	scalecolor(ar.rcol, zd);
150	>	zd *= (double)n;
151	>	scalecolor(ap->v, zd);
152	>	addcolor(ap->v, ar.rcol);
153	>	}
154	>	addcolor(hp->acol, ap->v);      /* add to our sum */
155	>	return(1);
156		}
157
158
159	+	/* Estimate errors based on ambient division differences */
160	+	static float *
161	+	getambdiffs(AMBHEMI *hp)
162	+	{
163	+	float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
164	+	float   *ep;
165	+	AMBSAMP *ap;
166	+	double  b, d2;
167	+	int     i, j;
168	+
169	+	if (earr == NULL)               /* out of memory? */
170	+	return(NULL);
171	+	/* compute squared neighbor diffs */
172	+	for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
173	+	for (j = 0; j < hp->ns; j++, ap++, ep++) {
174	+	b = bright(ap[0].v);
175	+	if (i) {                /* from above */
176	+	d2 = b - bright(ap[-hp->ns].v);
177	+	d2 *= d2;
178	+	ep[0] += d2;
179	+	ep[-hp->ns] += d2;
180	+	}
181	+	if (!j) continue;
182	+	/* from behind */
183	+	d2 = b - bright(ap[-1].v);
184	+	d2 *= d2;
185	+	ep[0] += d2;
186	+	ep[-1] += d2;
187	+	if (!i) continue;
188	+	/* diagonal */
189	+	d2 = b - bright(ap[-hp->ns-1].v);
190	+	d2 *= d2;
191	+	ep[0] += d2;
192	+	ep[-hp->ns-1] += d2;
193	+	}
194	+	/* correct for number of neighbors */
195	+	earr[0] *= 8./3.;
196	+	earr[hp->ns-1] *= 8./3.;
197	+	earr[(hp->ns-1)*hp->ns] *= 8./3.;
198	+	earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 8./3.;
199	+	for (i = 1; i < hp->ns-1; i++) {
200	+	earr[i*hp->ns] *= 8./5.;
201	+	earr[i*hp->ns + hp->ns-1] *= 8./5.;
202	+	}
203	+	for (j = 1; j < hp->ns-1; j++) {
204	+	earr[j] *= 8./5.;
205	+	earr[(hp->ns-1)*hp->ns + j] *= 8./5.;
206	+	}
207	+	return(earr);
208	+	}
209	+
210	+
211	+	/* Perform super-sampling on hemisphere (introduces bias) */
212	+	static void
213	+	ambsupersamp(AMBHEMI *hp, int cnt)
214	+	{
215	+	float   *earr = getambdiffs(hp);
216	+	double  e2rem = 0;
217	+	AMBSAMP *ap;
218	+	float   *ep;
219	+	int     i, j, n, nss;
220	+
221	+	if (earr == NULL)               /* just skip calc. if no memory */
222	+	return;
223	+	/* accumulate estimated variances */
224	+	for (ep = earr + hp->ns*hp->ns; ep > earr; )
225	+	e2rem += *--ep;
226	+	ep = earr;                      /* perform super-sampling */
227	+	for (ap = hp->sa, i = 0; i < hp->ns; i++)
228	+	for (j = 0; j < hp->ns; j++, ap++) {
229	+	if (e2rem <= FTINY)
230	+	goto done;      /* nothing left to do */
231	+	nss = *ep/e2rem*cnt + frandom();
232	+	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
233	+	--cnt;
234	+	e2rem -= *ep++;         /* update remainder */
235	+	}
236	+	done:
237	+	free(earr);
238	+	}
239	+
240	+
241	+	static AMBHEMI *
242	+	samp_hemi(                              /* sample indirect hemisphere */
243	+	COLOR   rcol,
244	+	RAY     *r,
245	+	double  wt
246	+	)
247	+	{
248	+	AMBHEMI *hp;
249	+	double  d;
250	+	int     n, i, j;
251	+	/* set number of divisions */
252	+	if (ambacc <= FTINY &&
253	+	wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight)))
254	+	wt = d;                 /* avoid ray termination */
255	+	n = sqrt(ambdiv * wt) + 0.5;
256	+	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
257	+	if (n < i)
258	+	n = i;
259	+	/* allocate sampling array */
260	+	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
261	+	if (hp == NULL)
262	+	error(SYSTEM, "out of memory in samp_hemi");
263	+	hp->rp = r;
264	+	hp->ns = n;
265	+	hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0;
266	+	memset(hp->sa, 0, sizeof(AMBSAMP)*n*n);
267	+	hp->sampOK = 0;
268	+	/* assign coefficient */
269	+	copycolor(hp->acoef, rcol);
270	+	d = 1.0/(n*n);
271	+	scalecolor(hp->acoef, d);
272	+	/* make tangent plane axes */
273	+	if (!getperpendicular(hp->ux, r->ron, 1))
274	+	error(CONSISTENCY, "bad ray direction in samp_hemi");
275	+	VCROSS(hp->uy, r->ron, hp->ux);
276	+	/* sample divisions */
277	+	for (i = hp->ns; i--; )
278	+	for (j = hp->ns; j--; )
279	+	hp->sampOK += ambsample(hp, i, j, 0);
280	+	copycolor(rcol, hp->acol);
281	+	if (!hp->sampOK) {              /* utter failure? */
282	+	free(hp);
283	+	return(NULL);
284	+	}
285	+	if (hp->sampOK < hp->ns*hp->ns) {
286	+	hp->sampOK *= -1;       /* soft failure */
287	+	return(hp);
288	+	}
289	+	n = ambssamp*wt + 0.5;
290	+	if (n > 8) {                    /* perform super-sampling? */
291	+	ambsupersamp(hp, n);
292	+	copycolor(rcol, hp->acol);
293	+	}
294	+	return(hp);                     /* all is well */
295	+	}
296	+
297	+
298	+	/* Return brightness of farthest ambient sample */
299	+	static double
300	+	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
301	+	{
302	+	if (hp->sa[n1].d <= hp->sa[n2].d) {
303	+	if (hp->sa[n1].d <= hp->sa[n3].d)
304	+	return(colval(hp->sa[n1].v,CIEY));
305	+	return(colval(hp->sa[n3].v,CIEY));
306	+	}
307	+	if (hp->sa[n2].d <= hp->sa[n3].d)
308	+	return(colval(hp->sa[n2].v,CIEY));
309	+	return(colval(hp->sa[n3].v,CIEY));
310	+	}
311	+
312	+
313		/* Compute vectors and coefficients for Hessian/gradient calcs */
314		static void
315	<	comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop)
315	>	comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1)
316		{
317	<	FVECT   vcp;
318	<	double  dot_e, dot_er, rdot_r, rdot_r1, J2;
145	<	int     i;
317	>	double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
318	>	int     ii;
319
320	<	VSUB(ftp->r_i, ap0, rop);
321	<	VSUB(ftp->r_i1, ap1, rop);
322	<	VSUB(ftp->e_i, ap1, ap0);
323	<	VCROSS(vcp, ftp->e_i, ftp->r_i);
324	<	ftp->nf = 1.0/DOT(vcp,vcp);
320	>	VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop);
321	>	VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop);
322	>	VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p);
323	>	VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
324	>	rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
325		dot_e = DOT(ftp->e_i,ftp->e_i);
326		dot_er = DOT(ftp->e_i, ftp->r_i);
327		rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i);
328		rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
329		ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) *
330	<	sqrt( ftp->nf );
330	>	sqrt( rdot_cp );
331		ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
332	<	dot_e*ftp->I1 )*0.5*ftp->nf;
332	>	dot_e*ftp->I1 )*0.5*rdot_cp;
333		J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
334	<	for (i = 3; i--; )
335	<	ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i];
334	>	for (ii = 3; ii--; )
335	>	ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
336		}
337
338
#	Line 180 | Line 353 | compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
353		static void
354		comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
355		{
356	<	FVECT   vcp;
356	>	FVECT   ncp;
357		FVECT   m1[3], m2[3], m3[3], m4[3];
358		double  d1, d2, d3, d4;
359		double  I3, J3, K3;
#	Line 190 | Line 363 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
363		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
364		d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
365		d4 = DOT(ftp->e_i, ftp->r_i);
366	<	I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 +
367	<	3.0/d3*ftp->I2 );
366	>	I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 )
367	>	/ ( 4.0*DOT(ftp->rcp,ftp->rcp) );
368		J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3;
369		K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3);
370		/* intermediate matrices */
371	<	VCROSS(vcp, nrm, ftp->e_i);
372	<	compose_matrix(m1, vcp, ftp->rI2_eJ2);
371	>	VCROSS(ncp, nrm, ftp->e_i);
372	>	compose_matrix(m1, ncp, ftp->rI2_eJ2);
373		compose_matrix(m2, ftp->r_i, ftp->r_i);
374		compose_matrix(m3, ftp->e_i, ftp->e_i);
375		compose_matrix(m4, ftp->r_i, ftp->e_i);
376	<	VCROSS(vcp, ftp->r_i, ftp->e_i);
204	<	d1 = DOT(nrm, vcp);
376	>	d1 = DOT(nrm, ftp->rcp);
377		d2 = -d1*ftp->I2;
378		d1 *= 2.0;
379		for (i = 3; i--; )              /* final matrix sum */
#	Line 209 | Line 381 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
381		hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] +
382		2.0*J3*m4[i][j] );
383		hess[i][j] += d2*(i==j);
384	<	hess[i][j] *= 1.0/PI;
384	>	hess[i][j] *= -1.0/PI;
385		}
386		}
387
#	Line 231 | Line 403 | rev_hessian(FVECT hess[3])
403		/* Add to radiometric Hessian from the given triangle */
404		static void
405		add2hessian(FVECT hess[3], FVECT ehess1[3],
406	<	FVECT ehess2[3], FVECT ehess3[3], COLORV v)
406	>	FVECT ehess2[3], FVECT ehess3[3], double v)
407		{
408		int     i, j;
409
#	Line 245 | Line 417 | add2hessian(FVECT hess[3], FVECT ehess1[3],
417		static void
418		comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
419		{
420	<	FVECT   vcp;
420	>	FVECT   ncp;
421		double  f1;
422		int     i;
423
424	<	VCROSS(vcp, ftp->r_i, ftp->r_i1);
425	<	f1 = 2.0*DOT(nrm, vcp);
254	<	VCROSS(vcp, nrm, ftp->e_i);
424	>	f1 = 2.0*DOT(nrm, ftp->rcp);
425	>	VCROSS(ncp, nrm, ftp->e_i);
426		for (i = 3; i--; )
427	<	grad[i] = (-0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] );
427	>	grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
428		}
429
430
#	Line 269 | Line 440 | rev_gradient(FVECT grad)
440
441		/* Add to displacement gradient from the given triangle */
442		static void
443	<	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
443	>	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
444		{
445		int     i;
446
#	Line 278 | Line 449 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
449		}
450
451
281	–	/* Return brightness of furthest ambient sample */
282	–	static COLORV
283	–	back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2,
284	–	struct s_ambsamp *ap3, FVECT orig)
285	–	{
286	–	COLORV  vback;
287	–	FVECT   vec;
288	–	double  d2, d2best;
289	–
290	–	VSUB(vec, ap1->p, orig);
291	–	d2best = DOT(vec,vec);
292	–	vback = colval(ap1->v,CIEY);
293	–	VSUB(vec, ap2->p, orig);
294	–	d2 = DOT(vec,vec);
295	–	if (d2 > d2best) {
296	–	d2best = d2;
297	–	vback = colval(ap2->v,CIEY);
298	–	}
299	–	VSUB(vec, ap3->p, orig);
300	–	d2 = DOT(vec,vec);
301	–	if (d2 > d2best)
302	–	return(colval(ap3->v,CIEY));
303	–	return(vback);
304	–	}
305	–
306	–
452		/* Compute anisotropic radii and eigenvector directions */
453	<	static int
453	>	static void
454		eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
455		{
456		double  hess2[2][2];
#	Line 321 | Line 466 | eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
466		hess2[0][1] = DOT(uv[0], b);
467		hess2[1][0] = DOT(uv[1], a);
468		hess2[1][1] = DOT(uv[1], b);
469	<	/* compute eigenvalues */
470	<	if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
471	<	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 ||
472	<	(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY ||
473	<	(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY )
474	<	error(INTERNAL, "bad eigenvalue calculation");
475	<
469	>	/* compute eigenvalue(s) */
470	>	i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
471	>	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]);
472	>	if (i == 1)                     /* double-root (circle) */
473	>	evalue[1] = evalue[0];
474	>	if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) |
475	>	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) {
476	>	ra[0] = ra[1] = maxarad;
477	>	return;
478	>	}
479		if (evalue[0] > evalue[1]) {
480		ra[0] = sqrt(sqrt(4.0/evalue[0]));
481		ra[1] = sqrt(sqrt(4.0/evalue[1]));
#	Line 385 | Line 533 | ambHessian(                            /* anisotropic radii & pos. gradient */
533		}
534		/* compute first row of edges */
535		for (j = 0; j < hp->ns-1; j++) {
536	<	comp_fftri(&fftr, ambsamp(hp,0,j).p,
389	<	ambsamp(hp,0,j+1).p, hp->rp->rop);
536	>	comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1));
537		if (hessrow != NULL)
538		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
539		if (gradrow != NULL)
#	Line 396 | Line 543 | ambHessian(                            /* anisotropic radii & pos. gradient */
543		for (i = 0; i < hp->ns-1; i++) {
544		FVECT       hesscol[3];     /* compute first vertical edge */
545		FVECT       gradcol;
546	<	comp_fftri(&fftr, ambsamp(hp,i,0).p,
400	<	ambsamp(hp,i+1,0).p, hp->rp->rop);
546	>	comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0));
547		if (hessrow != NULL)
548		comp_hessian(hesscol, &fftr, hp->rp->ron);
549		if (gradrow != NULL)
#	Line 405 | Line 551 | ambHessian(                            /* anisotropic radii & pos. gradient */
551		for (j = 0; j < hp->ns-1; j++) {
552		FVECT   hessdia[3];     /* compute triangle contributions */
553		FVECT   graddia;
554	<	COLORV  backg;
555	<	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
556	<	&ambsamp(hp,i+1,j), hp->rp->rop);
554	>	double  backg;
555	>	backg = back_ambval(hp, AI(hp,i,j),
556	>	AI(hp,i,j+1), AI(hp,i+1,j));
557		/* diagonal (inner) edge */
558	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
413	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
558	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j));
559		if (hessrow != NULL) {
560		comp_hessian(hessdia, &fftr, hp->rp->ron);
561		rev_hessian(hesscol);
562		add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
563		}
564	<	if (gradient != NULL) {
564	>	if (gradrow != NULL) {
565		comp_gradient(graddia, &fftr, hp->rp->ron);
566		rev_gradient(gradcol);
567		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
568		}
569		/* initialize edge in next row */
570	<	comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
426	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
570	>	comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j));
571		if (hessrow != NULL)
572		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
573		if (gradrow != NULL)
574		comp_gradient(gradrow[j], &fftr, hp->rp->ron);
575		/* new column edge & paired triangle */
576	<	backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
577	<	&ambsamp(hp,i+1,j), hp->rp->rop);
578	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
435	<	hp->rp->rop);
576	>	backg = back_ambval(hp, AI(hp,i+1,j+1),
577	>	AI(hp,i+1,j), AI(hp,i,j+1));
578	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1));
579		if (hessrow != NULL) {
580		comp_hessian(hesscol, &fftr, hp->rp->ron);
581		rev_hessian(hessdia);
#	Line 466 | Line 609 | ambHessian(                            /* anisotropic radii & pos. gradient */
609		static void
610		ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
611		{
612	<	struct s_ambsamp        *ap;
613	<	double                  dgsum[2];
614	<	int                     n;
615	<	FVECT                   vd;
616	<	double                  gfact;
612	>	AMBSAMP *ap;
613	>	double  dgsum[2];
614	>	int     n;
615	>	FVECT   vd;
616	>	double  gfact;
617
618		dgsum[0] = dgsum[1] = 0.0;      /* sum values times -tan(theta) */
619		for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
#	Line 478 | Line 621 | ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
621		VSUB(vd, ap->p, hp->rp->rop);
622		/* brightness over cosine factor */
623		gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
624	<	/* -sine = -proj_radius/vd_length */
625	<	dgsum[0] += DOT(uv[1], vd) * gfact;
626	<	dgsum[1] -= DOT(uv[0], vd) * gfact;
624	>	/* sine = proj_radius/vd_length */
625	>	dgsum[0] -= DOT(uv[1], vd) * gfact;
626	>	dgsum[1] += DOT(uv[0], vd) * gfact;
627		}
628		dg[0] = dgsum[0] / (hp->ns*hp->ns);
629		dg[1] = dgsum[1] / (hp->ns*hp->ns);
630		}
631
632
633	+	/* Compute potential light leak direction flags for cache value */
634	+	static uint32
635	+	ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
636	+	{
637	+	const double    max_d = 1.0/(minarad*ambacc + 0.001);
638	+	const double    ang_res = 0.5*PI/hp->ns;
639	+	const double    ang_step = ang_res/((int)(16/PI*ang_res) + 1.01);
640	+	double          avg_d = 0;
641	+	uint32          flgs = 0;
642	+	FVECT           vec;
643	+	double          u, v;
644	+	double          ang, a1;
645	+	int             i, j;
646	+	/* don't bother for a few samples */
647	+	if (hp->ns < 8)
648	+	return(0);
649	+	/* check distances overhead */
650	+	for (i = hp->ns*3/4; i-- > hp->ns>>2; )
651	+	for (j = hp->ns*3/4; j-- > hp->ns>>2; )
652	+	avg_d += ambsam(hp,i,j).d;
653	+	avg_d *= 4.0/(hp->ns*hp->ns);
654	+	if (avg_d*r0 >= 1.0)            /* ceiling too low for corral? */
655	+	return(0);
656	+	if (avg_d >= max_d)             /* insurance */
657	+	return(0);
658	+	/* else circle around perimeter */
659	+	for (i = 0; i < hp->ns; i++)
660	+	for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) {
661	+	AMBSAMP *ap = &ambsam(hp,i,j);
662	+	if ((ap->d <= FTINY) | (ap->d >= max_d))
663	+	continue;       /* too far or too near */
664	+	VSUB(vec, ap->p, hp->rp->rop);
665	+	u = DOT(vec, uv[0]);
666	+	v = DOT(vec, uv[1]);
667	+	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v)
668	+	continue;       /* occluder outside ellipse */
669	+	ang = atan2a(v, u);     /* else set direction flags */
670	+	for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step)
671	+	flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0)));
672	+	}
673	+	/* add low-angle incident (< 20deg) */
674	+	if (fabs(hp->rp->rod) <= 0.342) {
675	+	u = -DOT(hp->rp->rdir, uv[0]);
676	+	v = -DOT(hp->rp->rdir, uv[1]);
677	+	if ((r0*r0*u*u + r1*r1*v*v) > hp->rp->rot*hp->rp->rot) {
678	+	ang = atan2a(v, u);
679	+	ang += 2.*PI*(ang < 0);
680	+	ang *= 16/PI;
681	+	if ((ang < .5) | (ang >= 31.5))
682	+	flgs |= 0x80000001;
683	+	else
684	+	flgs |= 3L<<(int)(ang-.5);
685	+	}
686	+	}
687	+	return(flgs);
688	+	}
689	+
690	+
691		int
692		doambient(                              /* compute ambient component */
693		COLOR   rcol,                   /* input/output color */
#	Line 495 | Line 696 | doambient(                             /* compute ambient component */
696		FVECT   uv[2],                  /* returned (optional) */
697		float   ra[2],                  /* returned (optional) */
698		float   pg[2],                  /* returned (optional) */
699	<	float   dg[2]                   /* returned (optional) */
699	>	float   dg[2],                  /* returned (optional) */
700	>	uint32  *crlp                   /* returned (optional) */
701		)
702		{
703	<	AMBHEMI                 *hp = inithemi(rcol, r, wt);
704	<	int                     cnt = 0;
705	<	FVECT                   my_uv[2];
706	<	double                  d, acol[3];
707	<	struct s_ambsamp        *ap;
708	<	int                     i, j;
507	<	/* check/initialize */
508	<	if (hp == NULL)
509	<	return(0);
703	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
704	>	FVECT   my_uv[2];
705	>	double  d, K;
706	>	AMBSAMP *ap;
707	>	int     i;
708	>	/* clear return values */
709		if (uv != NULL)
710		memset(uv, 0, sizeof(FVECT)*2);
711		if (ra != NULL)
#	Line 515 | Line 714 | doambient(                             /* compute ambient component */
714		pg[0] = pg[1] = 0.0;
715		if (dg != NULL)
716		dg[0] = dg[1] = 0.0;
717	<	/* sample the hemisphere */
718	<	acol[0] = acol[1] = acol[2] = 0.0;
719	<	for (i = hp->ns; i--; )
720	<	for (j = hp->ns; j--; )
721	<	if ((ap = ambsample(hp, i, j)) != NULL) {
722	<	addcolor(acol, ap->v);
723	<	++cnt;
724	<	}
725	<	if (!cnt) {
527	<	setcolor(rcol, 0.0, 0.0, 0.0);
528	<	free(hp);
529	<	return(0);              /* no valid samples */
717	>	if (crlp != NULL)
718	>	*crlp = 0;
719	>	if (hp == NULL)                 /* sampling falure? */
720	>	return(0);
721	>
722	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) ||
723	>	(hp->sampOK < 0) | (hp->ns < 6)) {
724	>	free(hp);               /* Hessian not requested/possible */
725	>	return(-1);             /* value-only return value */
726		}
727	<	copycolor(rcol, acol);          /* final indirect irradiance/PI */
728	<	if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
729	<	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
730	<	free(hp);
731	<	return(-1);             /* no radius or gradient calc. */
727	>	if ((d = bright(rcol)) > FTINY) {       /* normalize Y values */
728	>	d = 0.99*(hp->ns*hp->ns)/d;
729	>	K = 0.01;
730	>	} else {                        /* or fall back on geometric Hessian */
731	>	K = 1.0;
732	>	pg = NULL;
733	>	dg = NULL;
734	>	crlp = NULL;
735		}
537	–	if (bright(acol) > FTINY)       /* normalize Y values */
538	–	d = cnt/bright(acol);
539	–	else
540	–	d = 0.0;
736		ap = hp->sa;                    /* relative Y channel from here on... */
737		for (i = hp->ns*hp->ns; i--; ap++)
738	<	colval(ap->v,CIEY) = bright(ap->v)*d + 0.01;
738	>	colval(ap->v,CIEY) = bright(ap->v)*d + K;
739
740		if (uv == NULL)                 /* make sure we have axis pointers */
741		uv = my_uv;
#	Line 551 | Line 746 | doambient(                             /* compute ambient component */
746		ambdirgrad(hp, uv, dg);
747
748		if (ra != NULL) {               /* scale/clamp radii */
749	+	if (pg != NULL) {
750	+	if (ra[0]*(d = fabs(pg[0])) > 1.0)
751	+	ra[0] = 1.0/d;
752	+	if (ra[1]*(d = fabs(pg[1])) > 1.0)
753	+	ra[1] = 1.0/d;
754	+	if (ra[0] > ra[1])
755	+	ra[0] = ra[1];
756	+	}
757		if (ra[0] < minarad) {
758		ra[0] = minarad;
759		if (ra[1] < minarad)
760		ra[1] = minarad;
761		}
762	<	ra[0] *= d = 1.0/sqrt(sqrt(wt));
762	>	ra[0] *= d = 1.0/sqrt(wt);
763		if ((ra[1] *= d) > 2.0*ra[0])
764		ra[1] = 2.0*ra[0];
765		if (ra[1] > maxarad) {
766		ra[1] = maxarad;
767		if (ra[0] > maxarad)
768		ra[0] = maxarad;
769	+	}
770	+	/* flag encroached directions */
771	+	if (crlp != NULL)
772	+	*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc);
773	+	if (pg != NULL) {       /* cap gradient if necessary */
774	+	d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];
775	+	if (d > 1.0) {
776	+	d = 1.0/sqrt(d);
777	+	pg[0] *= d;
778	+	pg[1] *= d;
779	+	}
780		}
781		}
782		free(hp);                       /* clean up and return */

Diff Legend

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