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.41 by greg, Wed Apr 30 23:38:58 2014 UTC vs.
Revision 2.49 by greg, Wed May 7 01:16:02 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 21 | Line 25 | static const char      RCSid[] = "$Id$";
25
26		extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
27
28	+	/* vertex direction bit positions */
29	+	#define VDB_xy  0
30	+	#define VDB_y   01
31	+	#define VDB_x   02
32	+	#define VDB_Xy  03
33	+	#define VDB_xY  04
34	+	#define VDB_X   05
35	+	#define VDB_Y   06
36	+	#define VDB_XY  07
37	+	/* get opposite vertex direction bit */
38	+	#define VDB_OPP(f)      (~(f) & 07)
39	+	/* adjacent triangle vertex flags */
40	+	static const int  adjacent_trifl[8] = {
41	+	0,                      /* forbidden diagonal */
42	+	1<<VDB_x|1<<VDB_y|1<<VDB_Xy,
43	+	1<<VDB_y|1<<VDB_x|1<<VDB_xY,
44	+	1<<VDB_y|1<<VDB_Xy|1<<VDB_X,
45	+	1<<VDB_x|1<<VDB_xY|1<<VDB_Y,
46	+	1<<VDB_Xy|1<<VDB_X|1<<VDB_Y,
47	+	1<<VDB_xY|1<<VDB_Y|1<<VDB_X,
48	+	0,                      /* forbidden diagonal */
49	+	};
50	+
51		typedef struct {
52	+	COLOR   v;              /* hemisphere sample value */
53	+	float   d;              /* reciprocal distance (1/rt) */
54	+	FVECT   p;              /* intersection point */
55	+	} AMBSAMP;              /* sample value */
56	+
57	+	typedef struct {
58		RAY     *rp;            /* originating ray sample */
59		FVECT   ux, uy;         /* tangent axis unit vectors */
60		int     ns;             /* number of samples per axis */
61		COLOR   acoef;          /* division contribution coefficient */
62	<	struct s_ambsamp {
30	<	COLOR   v;              /* hemisphere sample value */
31	<	FVECT   p;              /* intersection point */
32	<	} sa[1];                /* sample array (extends struct) */
62	>	AMBSAMP sa[1];          /* sample array (extends struct) */
63		}  AMBHEMI;             /* ambient sample hemisphere */
64
65	<	typedef struct s_ambsamp        AMBSAMP;
65	>	#define ambndx(h,i,j)   ((i)*(h)->ns + (j))
66	>	#define ambsam(h,i,j)   (h)->sa[ambndx(h,i,j)]
67
37	–	#define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
38	–
68		typedef struct {
69		FVECT   r_i, r_i1, e_i, rcp, rI2_eJ2;
70		double  I1, I2;
71	+	int     valid;
72		} FFTRI;                /* vectors and coefficients for Hessian calculation */
73
74
75	+	/* Get index for adjacent vertex */
76	+	static int
77	+	adjacent_verti(AMBHEMI *hp, int i, int j, int dbit)
78	+	{
79	+	int     i0 = i*hp->ns + j;
80	+
81	+	switch (dbit) {
82	+	case VDB_y:     return(i0 - hp->ns);
83	+	case VDB_x:     return(i0 - 1);
84	+	case VDB_Xy:    return(i0 - hp->ns + 1);
85	+	case VDB_xY:    return(i0 + hp->ns - 1);
86	+	case VDB_X:     return(i0 + 1);
87	+	case VDB_Y:     return(i0 + hp->ns);
88	+	/* the following should never occur */
89	+	case VDB_xy:    return(i0 - hp->ns - 1);
90	+	case VDB_XY:    return(i0 + hp->ns + 1);
91	+	}
92	+	return(-1);
93	+	}
94	+
95	+
96	+	/* Get vertex direction bit for the opposite edge to complete triangle */
97	+	static int
98	+	vdb_edge(int db1, int db2)
99	+	{
100	+	switch (db1) {
101	+	case VDB_x:     return(db2==VDB_y ? VDB_Xy : VDB_Y);
102	+	case VDB_y:     return(db2==VDB_x ? VDB_xY : VDB_X);
103	+	case VDB_X:     return(db2==VDB_Xy ? VDB_y : VDB_xY);
104	+	case VDB_Y:     return(db2==VDB_xY ? VDB_x : VDB_Xy);
105	+	case VDB_xY:    return(db2==VDB_x ? VDB_y : VDB_X);
106	+	case VDB_Xy:    return(db2==VDB_y ? VDB_x : VDB_Y);
107	+	}
108	+	error(INTERNAL, "forbidden diagonal in vdb_edge()");
109	+	return(-1);
110	+	}
111	+
112	+
113		static AMBHEMI *
114		inithemi(                       /* initialize sampling hemisphere */
115		COLOR   ac,
#	Line 88 | Line 156 | inithemi(                      /* initialize sampling hemisphere */
156		}
157
158
159	<	/* Prepare ambient division sample */
159	>	/* Sample ambient division and apply weighting coefficient */
160		static int
161	<	prepambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
161	>	getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
162		{
163		int     hlist[3], ii;
164		double  spt[2], zd;
#	Line 106 | Line 174 | prepambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n
174		scalecolor(arp->rcoef, 1./AVGREFL);
175		}
176		hlist[0] = hp->rp->rno;
177	<	hlist[1] = i;
178	<	hlist[2] = j;
177	>	hlist[1] = j;
178	>	hlist[2] = i;
179		multisamp(spt, 2, urand(ilhash(hlist,3)+n));
180		if (!n) {                       /* avoid border samples for n==0 */
181	<	if ((spt[0] < 0.1) | (spt[0] > 0.9))
181	>	if ((spt[0] < 0.1) | (spt[0] >= 0.9))
182		spt[0] = 0.1 + 0.8*frandom();
183	<	if ((spt[1] < 0.1) | (spt[1] > 0.9))
183	>	if ((spt[1] < 0.1) | (spt[1] >= 0.9))
184		spt[1] = 0.1 + 0.8*frandom();
185		}
186	<	SDsquare2disk(spt, (i+spt[0])/hp->ns, (j+spt[1])/hp->ns);
186	>	SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
187		zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
188		for (ii = 3; ii--; )
189		arp->rdir[ii] = spt[0]*hp->ux[ii] +
190		spt[1]*hp->uy[ii] +
191		zd*hp->rp->ron[ii];
192		checknorm(arp->rdir);
193	+	dimlist[ndims++] = ambndx(hp,i,j) + 90171;
194	+	rayvalue(arp);                  /* evaluate ray */
195	+	ndims--;                        /* apply coefficient */
196	+	multcolor(arp->rcol, arp->rcoef);
197		return(1);
198		}
199
200
201		static AMBSAMP *
202	<	ambsample(                              /* sample an ambient direction */
202	>	ambsample(                              /* initial ambient division sample */
203		AMBHEMI *hp,
204		int     i,
205		int     j
206		)
207		{
208	<	AMBSAMP *ap = &ambsamp(hp,i,j);
208	>	AMBSAMP *ap = &ambsam(hp,i,j);
209		RAY     ar;
210		/* generate hemispherical sample */
211	<	if (!prepambsamp(&ar, hp, i, j, 0))
212	<	goto badsample;
213	<	dimlist[ndims++] = i*hp->ns + j + 90171;
214	<	rayvalue(&ar);                  /* evaluate ray */
215	<	ndims--;
144	<	/* limit vertex distance */
211	>	if (!getambsamp(&ar, hp, i, j, 0) || ar.rt <= FTINY) {
212	>	memset(ap, 0, sizeof(AMBSAMP));
213	>	return(NULL);
214	>	}
215	>	ap->d = 1.0/ar.rt;              /* limit vertex distance */
216		if (ar.rt > 10.0*thescene.cusize)
217		ar.rt = 10.0*thescene.cusize;
147	–	else if (ar.rt <= FTINY)        /* should never happen! */
148	–	goto badsample;
218		VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
150	–	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
219		copycolor(ap->v, ar.rcol);
220		return(ap);
153	–	badsample:
154	–	setcolor(ap->v, 0., 0., 0.);
155	–	VCOPY(ap->p, hp->rp->rop);
156	–	return(NULL);
221		}
222
223
#	Line 161 | Line 225 | badsample:
225		static float *
226		getambdiffs(AMBHEMI *hp)
227		{
228	<	float   *earr = calloc(hp->ns*hp->ns, sizeof(float));
228	>	float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
229		float   *ep;
230	+	AMBSAMP *ap;
231		double  b, d2;
232		int     i, j;
233
234		if (earr == NULL)               /* out of memory? */
235		return(NULL);
236		/* compute squared neighbor diffs */
237	<	for (ep = earr, i = 0; i < hp->ns; i++)
238	<	for (j = 0; j < hp->ns; j++, ep++) {
239	<	b = bright(ambsamp(hp,i,j).v);
237	>	for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
238	>	for (j = 0; j < hp->ns; j++, ap++, ep++) {
239	>	b = bright(ap[0].v);
240		if (i) {                /* from above */
241	<	d2 = b - bright(ambsamp(hp,i-1,j).v);
241	>	d2 = b - bright(ap[-hp->ns].v);
242		d2 *= d2;
243		ep[0] += d2;
244		ep[-hp->ns] += d2;
245		}
246		if (j) {                /* from behind */
247	<	d2 = b - bright(ambsamp(hp,i,j-1).v);
247	>	d2 = b - bright(ap[-1].v);
248		d2 *= d2;
249		ep[0] += d2;
250		ep[-1] += d2;
#	Line 202 | Line 267 | getambdiffs(AMBHEMI *hp)
267		}
268
269
270	<	/* Perform super-sampling on hemisphere */
270	>	/* Perform super-sampling on hemisphere (introduces bias) */
271		static void
272		ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
273		{
#	Line 225 | Line 290 | ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
290		int     nss = *ep/e2sum*cnt + frandom();
291		setcolor(asum, 0., 0., 0.);
292		for (n = 1; n <= nss; n++) {
293	<	if (!prepambsamp(&ar, hp, i, j, n)) {
293	>	if (!getambsamp(&ar, hp, i, j, n)) {
294		nss = n-1;
295		break;
296		}
232	–	dimlist[ndims++] = i*hp->ns + j + 90171;
233	–	rayvalue(&ar);  /* evaluate super-sample */
234	–	ndims--;
235	–	multcolor(ar.rcol, ar.rcoef);
297		addcolor(asum, ar.rcol);
298		}
299		if (nss) {              /* update returned ambient value */
#	Line 248 | Line 309 | ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
309		}
310
311
312	+	/* Compute vertex flags, indicating farthest in each direction */
313	+	static uby8 *
314	+	vertex_flags(AMBHEMI *hp)
315	+	{
316	+	uby8    *vflags = (uby8 *)calloc(hp->ns*hp->ns, sizeof(uby8));
317	+	uby8    *vf;
318	+	AMBSAMP *ap;
319	+	int     i, j;
320	+
321	+	if (vflags == NULL)
322	+	error(SYSTEM, "out of memory in vertex_flags()");
323	+	vf = vflags;
324	+	ap = hp->sa;            /* compute farthest along first row */
325	+	for (j = 0; j < hp->ns-1; j++, vf++, ap++)
326	+	if (ap[0].d <= ap[1].d)
327	+	vf[0] |= 1<<VDB_X;
328	+	else
329	+	vf[1] |= 1<<VDB_x;
330	+	++vf; ++ap;
331	+	/* flag subsequent rows */
332	+	for (i = 1; i < hp->ns; i++) {
333	+	for (j = 0; j < hp->ns-1; j++, vf++, ap++) {
334	+	if (ap[0].d <= ap[-hp->ns].d)   /* row before */
335	+	vf[0] |= 1<<VDB_y;
336	+	else
337	+	vf[-hp->ns] |= 1<<VDB_Y;
338	+	if (ap[0].d <= ap[1-hp->ns].d)  /* diagonal we care about */
339	+	vf[0] |= 1<<VDB_Xy;
340	+	else
341	+	vf[1-hp->ns] |= 1<<VDB_xY;
342	+	if (ap[0].d <= ap[1].d)         /* column after */
343	+	vf[0] |= 1<<VDB_X;
344	+	else
345	+	vf[1] |= 1<<VDB_x;
346	+	}
347	+	if (ap[0].d <= ap[-hp->ns].d)       /* final column edge */
348	+	vf[0] |= 1<<VDB_y;
349	+	else
350	+	vf[-hp->ns] |= 1<<VDB_Y;
351	+	++vf; ++ap;
352	+	}
353	+	return(vflags);
354	+	}
355	+
356	+
357	+	/* Return brightness of farthest ambient sample */
358	+	static double
359	+	back_ambval(AMBHEMI *hp, int i, int j, int dbit1, int dbit2, const uby8 *vflags)
360	+	{
361	+	const int       v0 = ambndx(hp,i,j);
362	+	const int       tflags = (1<<dbit1 | 1<<dbit2);
363	+	int             v1, v2;
364	+
365	+	if ((vflags[v0] & tflags) == tflags)    /* is v0 the farthest? */
366	+	return(colval(hp->sa[v0].v,CIEY));
367	+	v1 = adjacent_verti(hp, i, j, dbit1);
368	+	if (vflags[v0] & 1<<dbit2)              /* v1 farthest if v0>v2 */
369	+	return(colval(hp->sa[v1].v,CIEY));
370	+	v2 = adjacent_verti(hp, i, j, dbit2);
371	+	if (vflags[v0] & 1<<dbit1)              /* v2 farthest if v0>v1 */
372	+	return(colval(hp->sa[v2].v,CIEY));
373	+	/* else check if v1>v2 */
374	+	if (vflags[v1] & 1<<vdb_edge(dbit1,dbit2))
375	+	return(colval(hp->sa[v1].v,CIEY));
376	+	return(colval(hp->sa[v2].v,CIEY));
377	+	}
378	+
379	+
380		/* Compute vectors and coefficients for Hessian/gradient calcs */
381		static void
382	<	comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop)
382	>	comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags)
383		{
384	<	double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
385	<	int     i;
384	>	const int       i0 = ambndx(hp,i,j);
385	>	double          rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
386	>	int             i1, ii;
387
388	<	VSUB(ftp->r_i, ap0, rop);
389	<	VSUB(ftp->r_i1, ap1, rop);
390	<	VSUB(ftp->e_i, ap1, ap0);
388	>	ftp->valid = 0;                 /* check if we can skip this edge */
389	>	ii = adjacent_trifl[dbit];
390	>	if ((vflags[i0] & ii) == ii)    /* cancels if vertex used as value */
391	>	return;
392	>	i1 = adjacent_verti(hp, i, j, dbit);
393	>	ii = adjacent_trifl[VDB_OPP(dbit)];
394	>	if ((vflags[i1] & ii) == ii)    /* on either end (for both triangles) */
395	>	return;
396	>	/* else go ahead with calculation */
397	>	VSUB(ftp->r_i, hp->sa[i0].p, hp->rp->rop);
398	>	VSUB(ftp->r_i1, hp->sa[i1].p, hp->rp->rop);
399	>	VSUB(ftp->e_i, hp->sa[i1].p, hp->sa[i0].p);
400		VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
401		rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
402		dot_e = DOT(ftp->e_i,ftp->e_i);
#	Line 269 | Line 408 | comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop
408		ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
409		dot_e*ftp->I1 )*0.5*rdot_cp;
410		J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
411	<	for (i = 3; i--; )
412	<	ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i];
411	>	for (ii = 3; ii--; )
412	>	ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
413	>	ftp->valid++;
414		}
415
416
#	Line 296 | Line 436 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
436		double  d1, d2, d3, d4;
437		double  I3, J3, K3;
438		int     i, j;
439	+
440	+	if (!ftp->valid) {              /* preemptive test */
441	+	memset(hess, 0, sizeof(FVECT)*3);
442	+	return;
443	+	}
444		/* compute intermediate coefficients */
445		d1 = 1.0/DOT(ftp->r_i,ftp->r_i);
446		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
#	Line 319 | Line 464 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
464		hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] +
465		2.0*J3*m4[i][j] );
466		hess[i][j] += d2*(i==j);
467	<	hess[i][j] *= 1.0/PI;
467	>	hess[i][j] *= -1.0/PI;
468		}
469		}
470
#	Line 341 | Line 486 | rev_hessian(FVECT hess[3])
486		/* Add to radiometric Hessian from the given triangle */
487		static void
488		add2hessian(FVECT hess[3], FVECT ehess1[3],
489	<	FVECT ehess2[3], FVECT ehess3[3], COLORV v)
489	>	FVECT ehess2[3], FVECT ehess3[3], double v)
490		{
491		int     i, j;
492
#	Line 359 | Line 504 | comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
504		double  f1;
505		int     i;
506
507	+	if (!ftp->valid) {              /* preemptive test */
508	+	memset(grad, 0, sizeof(FVECT));
509	+	return;
510	+	}
511		f1 = 2.0*DOT(nrm, ftp->rcp);
512		VCROSS(ncp, nrm, ftp->e_i);
513		for (i = 3; i--; )
514	<	grad[i] = (-0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
514	>	grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
515		}
516
517
#	Line 378 | Line 527 | rev_gradient(FVECT grad)
527
528		/* Add to displacement gradient from the given triangle */
529		static void
530	<	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
530	>	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
531		{
532		int     i;
533
#	Line 387 | Line 536 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
536		}
537
538
390	–	/* Return brightness of furthest ambient sample */
391	–	static COLORV
392	–	back_ambval(AMBSAMP *ap1, AMBSAMP *ap2, AMBSAMP *ap3, FVECT orig)
393	–	{
394	–	COLORV  vback;
395	–	FVECT   vec;
396	–	double  d2, d2best;
397	–
398	–	VSUB(vec, ap1->p, orig);
399	–	d2best = DOT(vec,vec);
400	–	vback = colval(ap1->v,CIEY);
401	–	VSUB(vec, ap2->p, orig);
402	–	d2 = DOT(vec,vec);
403	–	if (d2 > d2best) {
404	–	d2best = d2;
405	–	vback = colval(ap2->v,CIEY);
406	–	}
407	–	VSUB(vec, ap3->p, orig);
408	–	d2 = DOT(vec,vec);
409	–	if (d2 > d2best)
410	–	return(colval(ap3->v,CIEY));
411	–	return(vback);
412	–	}
413	–
414	–
539		/* Compute anisotropic radii and eigenvector directions */
540		static int
541		eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
#	Line 472 | Line 596 | ambHessian(                            /* anisotropic radii & pos. gradient */
596		static char     memerrmsg[] = "out of memory in ambHessian()";
597		FVECT           (*hessrow)[3] = NULL;
598		FVECT           *gradrow = NULL;
599	+	uby8            *vflags;
600		FVECT           hessian[3];
601		FVECT           gradient;
602		FFTRI           fftr;
#	Line 493 | Line 618 | ambHessian(                            /* anisotropic radii & pos. gradient */
618		error(SYSTEM, memerrmsg);
619		memset(gradient, 0, sizeof(gradient));
620		}
621	+	/* get vertex position flags */
622	+	vflags = vertex_flags(hp);
623		/* compute first row of edges */
624		for (j = 0; j < hp->ns-1; j++) {
625	<	comp_fftri(&fftr, ambsamp(hp,0,j).p,
499	<	ambsamp(hp,0,j+1).p, hp->rp->rop);
625	>	comp_fftri(&fftr, hp, 0, j, VDB_X, vflags);
626		if (hessrow != NULL)
627		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
628		if (gradrow != NULL)
#	Line 506 | Line 632 | ambHessian(                            /* anisotropic radii & pos. gradient */
632		for (i = 0; i < hp->ns-1; i++) {
633		FVECT       hesscol[3];     /* compute first vertical edge */
634		FVECT       gradcol;
635	<	comp_fftri(&fftr, ambsamp(hp,i,0).p,
510	<	ambsamp(hp,i+1,0).p, hp->rp->rop);
635	>	comp_fftri(&fftr, hp, i, 0, VDB_Y, vflags);
636		if (hessrow != NULL)
637		comp_hessian(hesscol, &fftr, hp->rp->ron);
638		if (gradrow != NULL)
#	Line 515 | Line 640 | ambHessian(                            /* anisotropic radii & pos. gradient */
640		for (j = 0; j < hp->ns-1; j++) {
641		FVECT   hessdia[3];     /* compute triangle contributions */
642		FVECT   graddia;
643	<	COLORV  backg;
644	<	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
520	<	&ambsamp(hp,i+1,j), hp->rp->rop);
643	>	double  backg;
644	>	backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags);
645		/* diagonal (inner) edge */
646	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
523	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
646	>	comp_fftri(&fftr, hp, i, j+1, VDB_xY, vflags);
647		if (hessrow != NULL) {
648		comp_hessian(hessdia, &fftr, hp->rp->ron);
649		rev_hessian(hesscol);
#	Line 532 | Line 655 | ambHessian(                            /* anisotropic radii & pos. gradient */
655		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
656		}
657		/* initialize edge in next row */
658	<	comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
536	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
658	>	comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags);
659		if (hessrow != NULL)
660		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
661		if (gradrow != NULL)
662		comp_gradient(gradrow[j], &fftr, hp->rp->ron);
663		/* new column edge & paired triangle */
664	<	backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
665	<	&ambsamp(hp,i+1,j), hp->rp->rop);
544	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
545	<	hp->rp->rop);
664	>	backg = back_ambval(hp, i+1, j+1, VDB_x, VDB_y, vflags);
665	>	comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags);
666		if (hessrow != NULL) {
667		comp_hessian(hesscol, &fftr, hp->rp->ron);
668		rev_hessian(hessdia);
#	Line 562 | Line 682 | ambHessian(                            /* anisotropic radii & pos. gradient */
682		/* release row buffers */
683		if (hessrow != NULL) free(hessrow);
684		if (gradrow != NULL) free(gradrow);
685	+	free(vflags);
686
687		if (ra != NULL)                 /* extract eigenvectors & radii */
688		eigenvectors(uv, ra, hessian);
#	Line 597 | Line 718 | ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
718		}
719
720
721	+	/* Compute potential light leak direction flags for cache value */
722	+	static uint32
723	+	ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
724	+	{
725	+	uint32  flgs = 0;
726	+	int     i, j;
727	+	/* circle around perimeter */
728	+	for (i = 0; i < hp->ns; i++)
729	+	for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) {
730	+	AMBSAMP *ap = &ambsam(hp,i,j);
731	+	FVECT   vec;
732	+	double  u, v;
733	+	double  ang;
734	+	int     abp;
735	+	if (ap->d <= FTINY)
736	+	continue;
737	+	VSUB(vec, ap->p, hp->rp->rop);
738	+	u = DOT(vec, uv[0]) * ap->d;
739	+	v = DOT(vec, uv[1]) * ap->d;
740	+	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0)
741	+	continue;       /* occluder outside ellipse */
742	+	ang = atan2a(v, u);     /* else set direction flags */
743	+	ang += 2.0*PI*(ang < 0);
744	+	ang *= 16./PI;
745	+	if ((ang < .5) | (ang >= 31.5))
746	+	flgs |= 0x80000001;
747	+	else
748	+	flgs |= 3L<<(int)(ang-.5);
749	+	}
750	+	return(flgs);
751	+	}
752	+
753	+
754		int
755		doambient(                              /* compute ambient component */
756		COLOR   rcol,                   /* input/output color */
#	Line 605 | Line 759 | doambient(                             /* compute ambient component */
759		FVECT   uv[2],                  /* returned (optional) */
760		float   ra[2],                  /* returned (optional) */
761		float   pg[2],                  /* returned (optional) */
762	<	float   dg[2]                   /* returned (optional) */
762	>	float   dg[2],                  /* returned (optional) */
763	>	uint32  *crlp                   /* returned (optional) */
764		)
765		{
766		AMBHEMI *hp = inithemi(rcol, r, wt);
767	<	int     cnt = 0;
767	>	int     cnt;
768		FVECT   my_uv[2];
769		double  d, K, acol[3];
770		AMBSAMP *ap;
#	Line 625 | Line 780 | doambient(                             /* compute ambient component */
780		pg[0] = pg[1] = 0.0;
781		if (dg != NULL)
782		dg[0] = dg[1] = 0.0;
783	+	if (crlp != NULL)
784	+	*crlp = 0;
785		/* sample the hemisphere */
786		acol[0] = acol[1] = acol[2] = 0.0;
787	+	cnt = 0;
788		for (i = hp->ns; i--; )
789		for (j = hp->ns; j--; )
790		if ((ap = ambsample(hp, i, j)) != NULL) {
#	Line 651 | Line 809 | doambient(                             /* compute ambient component */
809		free(hp);
810		return(-1);             /* no radius or gradient calc. */
811		}
812	<	if (bright(acol) > FTINY) {     /* normalize Y values */
813	<	d = 0.99*cnt/bright(acol);
812	>	if ((d = bright(acol)) > FTINY) {       /* normalize Y values */
813	>	d = 0.99*(hp->ns*hp->ns)/d;
814		K = 0.01;
815	<	} else {                        /* geometric Hessian fall-back */
658	<	d = 0.0;
815	>	} else {                        /* or fall back on geometric Hessian */
816		K = 1.0;
817		pg = NULL;
818		dg = NULL;
#	Line 694 | Line 851 | doambient(                             /* compute ambient component */
851		if (ra[0] > maxarad)
852		ra[0] = maxarad;
853		}
854	+	if (crlp != NULL)       /* flag encroached directions */
855	+	*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc);
856		if (pg != NULL) {       /* cap gradient if necessary */
857		d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];
858		if (d > 1.0) {

Diff Legend

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