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.20 by greg, Fri Oct 28 16:16:33 2005 UTC vs.
Revision 2.84 by greg, Tue Feb 26 00:37:54 2019 UTC

#	Line 4 | Line 4 | static const char      RCSid[] = "$Id$";
4		/*
5		* Routines to compute "ambient" values using Monte Carlo
6		*
7	+	*  Hessian calculations based on "Practical Hessian-Based Error Control
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
18		#include "copyright.h"
19
20		#include  "ray.h"
13	–
21		#include  "ambient.h"
15	–
22		#include  "random.h"
23
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 */
31	+	FVECT   p;              /* intersection point */
32	+	} AMBSAMP;              /* sample value */
33	+
34	+	typedef struct {
35	+	RAY     *rp;            /* originating ray sample */
36	+	int     ns;             /* number of samples per axis */
37	+	int     sampOK;         /* acquired full sample set? */
38	+	COLOR   acoef;          /* division contribution coefficient */
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 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, rcp, rI2_eJ2;
49	+	double  I1, I2;
50	+	} FFTRI;                /* vectors and coefficients for Hessian calculation */
51	+
52	+
53	+	static int
54	+	ambcollision(                           /* proposed direciton collides? */
55	+	AMBHEMI *hp,
56	+	int     i,
57	+	int     j,
58	+	FVECT   dv
59	+	)
60	+	{
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 int
91	+	ambsample(                              /* initial ambient division sample */
92	+	AMBHEMI *hp,
93	+	int     i,
94	+	int     j,
95	+	int     n
96	+	)
97	+	{
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	+	return(0);
110	+	if (ambacc > FTINY) {
111	+	multcolor(ar.rcoef, hp->acoef);
112	+	scalecolor(ar.rcoef, 1./AVGREFL);
113	+	}
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	+	/* 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	+	zd = raydistance(&ar);
135	+	if (zd <= FTINY)
136	+	return(0);              /* should never happen */
137	+	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
138	+	if (zd*ap->d < 1.0)             /* new/closer distance? */
139	+	ap->d = 1.0/zd;
140	+	if (!n) {                       /* record first vertex & value */
141	+	if (zd > 10.0*thescene.cusize + 1000.)
142	+	zd = 10.0*thescene.cusize + 1000.;
143	+	VSUM(ap->p, ar.rorg, ar.rdir, zd);
144	+	copycolor(ap->v, ar.rcol);
145	+	} else {                        /* else update recorded value */
146	+	hp->acol[RED] -= colval(ap->v,RED);
147	+	hp->acol[GRN] -= colval(ap->v,GRN);
148	+	hp->acol[BLU] -= colval(ap->v,BLU);
149	+	zd = 1.0/(double)(n+1);
150	+	scalecolor(ar.rcol, zd);
151	+	zd *= (double)n;
152	+	scalecolor(ap->v, zd);
153	+	addcolor(ap->v, ar.rcol);
154	+	}
155	+	addcolor(hp->acol, ap->v);      /* add to our sum */
156	+	return(1);
157	+	}
158	+
159	+
160	+	/* Estimate variance based on ambient division differences */
161	+	static float *
162	+	getambdiffs(AMBHEMI *hp)
163	+	{
164	+	const double    normf = 1./bright(hp->acoef);
165	+	float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
166	+	float   *ep;
167	+	AMBSAMP *ap;
168	+	double  b, b1, d2;
169	+	int     i, j;
170	+
171	+	if (earr == NULL)               /* out of memory? */
172	+	return(NULL);
173	+	/* sum squared neighbor diffs */
174	+	for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
175	+	for (j = 0; j < hp->ns; j++, ap++, ep++) {
176	+	b = bright(ap[0].v);
177	+	if (i) {                /* from above */
178	+	b1 = bright(ap[-hp->ns].v);
179	+	d2 = b - b1;
180	+	d2 *= d2*normf/(b + b1);
181	+	ep[0] += d2;
182	+	ep[-hp->ns] += d2;
183	+	}
184	+	if (!j) continue;
185	+	/* from behind */
186	+	b1 = bright(ap[-1].v);
187	+	d2 = b - b1;
188	+	d2 *= d2*normf/(b + b1);
189	+	ep[0] += d2;
190	+	ep[-1] += d2;
191	+	if (!i) continue;
192	+	/* diagonal */
193	+	b1 = bright(ap[-hp->ns-1].v);
194	+	d2 = b - b1;
195	+	d2 *= d2*normf/(b + b1);
196	+	ep[0] += d2;
197	+	ep[-hp->ns-1] += d2;
198	+	}
199	+	/* correct for number of neighbors */
200	+	earr[0] *= 8./3.;
201	+	earr[hp->ns-1] *= 8./3.;
202	+	earr[(hp->ns-1)*hp->ns] *= 8./3.;
203	+	earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 8./3.;
204	+	for (i = 1; i < hp->ns-1; i++) {
205	+	earr[i*hp->ns] *= 8./5.;
206	+	earr[i*hp->ns + hp->ns-1] *= 8./5.;
207	+	}
208	+	for (j = 1; j < hp->ns-1; j++) {
209	+	earr[j] *= 8./5.;
210	+	earr[(hp->ns-1)*hp->ns + j] *= 8./5.;
211	+	}
212	+	return(earr);
213	+	}
214	+
215	+
216	+	/* Perform super-sampling on hemisphere (introduces bias) */
217	+	static void
218	+	ambsupersamp(AMBHEMI *hp, int cnt)
219	+	{
220	+	float   *earr = getambdiffs(hp);
221	+	double  e2rem = 0;
222	+	float   *ep;
223	+	int     i, j, n, nss;
224	+
225	+	if (earr == NULL)               /* just skip calc. if no memory */
226	+	return;
227	+	/* accumulate estimated variances */
228	+	for (ep = earr + hp->ns*hp->ns; ep > earr; )
229	+	e2rem += *--ep;
230	+	ep = earr;                      /* perform super-sampling */
231	+	for (i = 0; i < hp->ns; i++)
232	+	for (j = 0; j < hp->ns; j++) {
233	+	if (e2rem <= FTINY)
234	+	goto done;      /* nothing left to do */
235	+	nss = *ep/e2rem*cnt + frandom();
236	+	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
237	+	if (!--cnt) goto done;
238	+	e2rem -= *ep++;         /* update remainder */
239	+	}
240	+	done:
241	+	free(earr);
242	+	}
243	+
244	+
245	+	static AMBHEMI *
246	+	samp_hemi(                              /* sample indirect hemisphere */
247	+	COLOR   rcol,
248	+	RAY     *r,
249	+	double  wt
250	+	)
251	+	{
252	+	AMBHEMI *hp;
253	+	double  d;
254	+	int     n, i, j;
255	+	/* insignificance check */
256	+	if (bright(rcol) <= FTINY)
257	+	return(NULL);
258	+	/* set number of divisions */
259	+	if (ambacc <= FTINY &&
260	+	wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight)))
261	+	wt = d;                 /* avoid ray termination */
262	+	n = sqrt(ambdiv * wt) + 0.5;
263	+	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
264	+	if (n < i)
265	+	n = i;
266	+	/* allocate sampling array */
267	+	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
268	+	if (hp == NULL)
269	+	error(SYSTEM, "out of memory in samp_hemi");
270	+	hp->rp = r;
271	+	hp->ns = n;
272	+	hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0;
273	+	memset(hp->sa, 0, sizeof(AMBSAMP)*n*n);
274	+	hp->sampOK = 0;
275	+	/* assign coefficient */
276	+	copycolor(hp->acoef, rcol);
277	+	d = 1.0/(n*n);
278	+	scalecolor(hp->acoef, d);
279	+	/* make tangent plane axes */
280	+	if (!getperpendicular(hp->ux, r->ron, 1))
281	+	error(CONSISTENCY, "bad ray direction in samp_hemi");
282	+	VCROSS(hp->uy, r->ron, hp->ux);
283	+	/* sample divisions */
284	+	for (i = hp->ns; i--; )
285	+	for (j = hp->ns; j--; )
286	+	hp->sampOK += ambsample(hp, i, j, 0);
287	+	copycolor(rcol, hp->acol);
288	+	if (!hp->sampOK) {              /* utter failure? */
289	+	free(hp);
290	+	return(NULL);
291	+	}
292	+	if (hp->sampOK < hp->ns*hp->ns) {
293	+	hp->sampOK *= -1;       /* soft failure */
294	+	return(hp);
295	+	}
296	+	if (hp->sampOK < 64)
297	+	return(hp);             /* insufficient for super-sampling */
298	+	n = ambssamp*wt + 0.5;
299	+	if (n > 8) {                    /* perform super-sampling? */
300	+	ambsupersamp(hp, n);
301	+	copycolor(rcol, hp->acol);
302	+	}
303	+	return(hp);                     /* all is well */
304	+	}
305	+
306	+
307	+	/* Return brightness of farthest ambient sample */
308	+	static double
309	+	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
310	+	{
311	+	if (hp->sa[n1].d <= hp->sa[n2].d) {
312	+	if (hp->sa[n1].d <= hp->sa[n3].d)
313	+	return(colval(hp->sa[n1].v,CIEY));
314	+	return(colval(hp->sa[n3].v,CIEY));
315	+	}
316	+	if (hp->sa[n2].d <= hp->sa[n3].d)
317	+	return(colval(hp->sa[n2].v,CIEY));
318	+	return(colval(hp->sa[n3].v,CIEY));
319	+	}
320	+
321	+
322	+	/* Compute vectors and coefficients for Hessian/gradient calcs */
323	+	static void
324	+	comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1)
325	+	{
326	+	double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
327	+	int     ii;
328	+
329	+	VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop);
330	+	VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop);
331	+	VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p);
332	+	VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
333	+	rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
334	+	dot_e = DOT(ftp->e_i,ftp->e_i);
335	+	dot_er = DOT(ftp->e_i, ftp->r_i);
336	+	rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i);
337	+	rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
338	+	ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) *
339	+	sqrt( rdot_cp );
340	+	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
341	+	dot_e*ftp->I1 )*0.5*rdot_cp;
342	+	J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
343	+	for (ii = 3; ii--; )
344	+	ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
345	+	}
346	+
347	+
348	+	/* Compose 3x3 matrix from two vectors */
349	+	static void
350	+	compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
351	+	{
352	+	mat[0][0] = 2.0*va[0]*vb[0];
353	+	mat[1][1] = 2.0*va[1]*vb[1];
354	+	mat[2][2] = 2.0*va[2]*vb[2];
355	+	mat[0][1] = mat[1][0] = va[0]*vb[1] + va[1]*vb[0];
356	+	mat[0][2] = mat[2][0] = va[0]*vb[2] + va[2]*vb[0];
357	+	mat[1][2] = mat[2][1] = va[1]*vb[2] + va[2]*vb[1];
358	+	}
359	+
360	+
361	+	/* Compute partial 3x3 Hessian matrix for edge */
362	+	static void
363	+	comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
364	+	{
365	+	FVECT   ncp;
366	+	FVECT   m1[3], m2[3], m3[3], m4[3];
367	+	double  d1, d2, d3, d4;
368	+	double  I3, J3, K3;
369	+	int     i, j;
370	+	/* compute intermediate coefficients */
371	+	d1 = 1.0/DOT(ftp->r_i,ftp->r_i);
372	+	d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
373	+	d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
374	+	d4 = DOT(ftp->e_i, ftp->r_i);
375	+	I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 )
376	+	/ ( 4.0*DOT(ftp->rcp,ftp->rcp) );
377	+	J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3;
378	+	K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3);
379	+	/* intermediate matrices */
380	+	VCROSS(ncp, nrm, ftp->e_i);
381	+	compose_matrix(m1, ncp, ftp->rI2_eJ2);
382	+	compose_matrix(m2, ftp->r_i, ftp->r_i);
383	+	compose_matrix(m3, ftp->e_i, ftp->e_i);
384	+	compose_matrix(m4, ftp->r_i, ftp->e_i);
385	+	d1 = DOT(nrm, ftp->rcp);
386	+	d2 = -d1*ftp->I2;
387	+	d1 *= 2.0;
388	+	for (i = 3; i--; )              /* final matrix sum */
389	+	for (j = 3; j--; ) {
390	+	hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] +
391	+	2.0*J3*m4[i][j] );
392	+	hess[i][j] += d2*(i==j);
393	+	hess[i][j] *= -1.0/PI;
394	+	}
395	+	}
396	+
397	+
398	+	/* Reverse hessian calculation result for edge in other direction */
399	+	static void
400	+	rev_hessian(FVECT hess[3])
401	+	{
402	+	int     i;
403	+
404	+	for (i = 3; i--; ) {
405	+	hess[i][0] = -hess[i][0];
406	+	hess[i][1] = -hess[i][1];
407	+	hess[i][2] = -hess[i][2];
408	+	}
409	+	}
410	+
411	+
412	+	/* Add to radiometric Hessian from the given triangle */
413	+	static void
414	+	add2hessian(FVECT hess[3], FVECT ehess1[3],
415	+	FVECT ehess2[3], FVECT ehess3[3], double v)
416	+	{
417	+	int     i, j;
418	+
419	+	for (i = 3; i--; )
420	+	for (j = 3; j--; )
421	+	hess[i][j] += v*( ehess1[i][j] + ehess2[i][j] + ehess3[i][j] );
422	+	}
423	+
424	+
425	+	/* Compute partial displacement form factor gradient for edge */
426	+	static void
427	+	comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
428	+	{
429	+	FVECT   ncp;
430	+	double  f1;
431	+	int     i;
432	+
433	+	f1 = 2.0*DOT(nrm, ftp->rcp);
434	+	VCROSS(ncp, nrm, ftp->e_i);
435	+	for (i = 3; i--; )
436	+	grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
437	+	}
438	+
439	+
440	+	/* Reverse gradient calculation result for edge in other direction */
441	+	static void
442	+	rev_gradient(FVECT grad)
443	+	{
444	+	grad[0] = -grad[0];
445	+	grad[1] = -grad[1];
446	+	grad[2] = -grad[2];
447	+	}
448	+
449	+
450	+	/* Add to displacement gradient from the given triangle */
451	+	static void
452	+	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
453	+	{
454	+	int     i;
455	+
456	+	for (i = 3; i--; )
457	+	grad[i] += v*( egrad1[i] + egrad2[i] + egrad3[i] );
458	+	}
459	+
460	+
461	+	/* Compute anisotropic radii and eigenvector directions */
462	+	static void
463	+	eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
464	+	{
465	+	double  hess2[2][2];
466	+	FVECT   a, b;
467	+	double  evalue[2], slope1, xmag1;
468	+	int     i;
469	+	/* project Hessian to sample plane */
470	+	for (i = 3; i--; ) {
471	+	a[i] = DOT(hessian[i], uv[0]);
472	+	b[i] = DOT(hessian[i], uv[1]);
473	+	}
474	+	hess2[0][0] = DOT(uv[0], a);
475	+	hess2[0][1] = DOT(uv[0], b);
476	+	hess2[1][0] = DOT(uv[1], a);
477	+	hess2[1][1] = DOT(uv[1], b);
478	+	/* compute eigenvalue(s) */
479	+	i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
480	+	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]);
481	+	if (i == 1)                     /* double-root (circle) */
482	+	evalue[1] = evalue[0];
483	+	if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) |
484	+	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) {
485	+	ra[0] = ra[1] = maxarad;
486	+	return;
487	+	}
488	+	if (evalue[0] > evalue[1]) {
489	+	ra[0] = sqrt(sqrt(4.0/evalue[0]));
490	+	ra[1] = sqrt(sqrt(4.0/evalue[1]));
491	+	slope1 = evalue[1];
492	+	} else {
493	+	ra[0] = sqrt(sqrt(4.0/evalue[1]));
494	+	ra[1] = sqrt(sqrt(4.0/evalue[0]));
495	+	slope1 = evalue[0];
496	+	}
497	+	/* compute unit eigenvectors */
498	+	if (fabs(hess2[0][1]) <= FTINY)
499	+	return;                 /* uv OK as is */
500	+	slope1 = (slope1 - hess2[0][0]) / hess2[0][1];
501	+	xmag1 = sqrt(1.0/(1.0 + slope1*slope1));
502	+	for (i = 3; i--; ) {
503	+	b[i] = xmag1*uv[0][i] + slope1*xmag1*uv[1][i];
504	+	a[i] = slope1*xmag1*uv[0][i] - xmag1*uv[1][i];
505	+	}
506	+	VCOPY(uv[0], a);
507	+	VCOPY(uv[1], b);
508	+	}
509	+
510	+
511	+	static void
512	+	ambHessian(                             /* anisotropic radii & pos. gradient */
513	+	AMBHEMI *hp,
514	+	FVECT   uv[2],                  /* returned */
515	+	float   ra[2],                  /* returned (optional) */
516	+	float   pg[2]                   /* returned (optional) */
517	+	)
518	+	{
519	+	static char     memerrmsg[] = "out of memory in ambHessian()";
520	+	FVECT           (*hessrow)[3] = NULL;
521	+	FVECT           *gradrow = NULL;
522	+	FVECT           hessian[3];
523	+	FVECT           gradient;
524	+	FFTRI           fftr;
525	+	int             i, j;
526	+	/* be sure to assign unit vectors */
527	+	VCOPY(uv[0], hp->ux);
528	+	VCOPY(uv[1], hp->uy);
529	+	/* clock-wise vertex traversal from sample POV */
530	+	if (ra != NULL) {               /* initialize Hessian row buffer */
531	+	hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*(hp->ns-1));
532	+	if (hessrow == NULL)
533	+	error(SYSTEM, memerrmsg);
534	+	memset(hessian, 0, sizeof(hessian));
535	+	} else if (pg == NULL)          /* bogus call? */
536	+	return;
537	+	if (pg != NULL) {               /* initialize form factor row buffer */
538	+	gradrow = (FVECT *)malloc(sizeof(FVECT)*(hp->ns-1));
539	+	if (gradrow == NULL)
540	+	error(SYSTEM, memerrmsg);
541	+	memset(gradient, 0, sizeof(gradient));
542	+	}
543	+	/* compute first row of edges */
544	+	for (j = 0; j < hp->ns-1; j++) {
545	+	comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1));
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	+	}
551	+	/* sum each row of triangles */
552	+	for (i = 0; i < hp->ns-1; i++) {
553	+	FVECT       hesscol[3];     /* compute first vertical edge */
554	+	FVECT       gradcol;
555	+	comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0));
556	+	if (hessrow != NULL)
557	+	comp_hessian(hesscol, &fftr, hp->rp->ron);
558	+	if (gradrow != NULL)
559	+	comp_gradient(gradcol, &fftr, hp->rp->ron);
560	+	for (j = 0; j < hp->ns-1; j++) {
561	+	FVECT   hessdia[3];     /* compute triangle contributions */
562	+	FVECT   graddia;
563	+	double  backg;
564	+	backg = back_ambval(hp, AI(hp,i,j),
565	+	AI(hp,i,j+1), AI(hp,i+1,j));
566	+	/* diagonal (inner) edge */
567	+	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j));
568	+	if (hessrow != NULL) {
569	+	comp_hessian(hessdia, &fftr, hp->rp->ron);
570	+	rev_hessian(hesscol);
571	+	add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
572	+	}
573	+	if (gradrow != NULL) {
574	+	comp_gradient(graddia, &fftr, hp->rp->ron);
575	+	rev_gradient(gradcol);
576	+	add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
577	+	}
578	+	/* initialize edge in next row */
579	+	comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j));
580	+	if (hessrow != NULL)
581	+	comp_hessian(hessrow[j], &fftr, hp->rp->ron);
582	+	if (gradrow != NULL)
583	+	comp_gradient(gradrow[j], &fftr, hp->rp->ron);
584	+	/* new column edge & paired triangle */
585	+	backg = back_ambval(hp, AI(hp,i+1,j+1),
586	+	AI(hp,i+1,j), AI(hp,i,j+1));
587	+	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1));
588	+	if (hessrow != NULL) {
589	+	comp_hessian(hesscol, &fftr, hp->rp->ron);
590	+	rev_hessian(hessdia);
591	+	add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
592	+	if (i < hp->ns-2)
593	+	rev_hessian(hessrow[j]);
594	+	}
595	+	if (gradrow != NULL) {
596	+	comp_gradient(gradcol, &fftr, hp->rp->ron);
597	+	rev_gradient(graddia);
598	+	add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
599	+	if (i < hp->ns-2)
600	+	rev_gradient(gradrow[j]);
601	+	}
602	+	}
603	+	}
604	+	/* release row buffers */
605	+	if (hessrow != NULL) free(hessrow);
606	+	if (gradrow != NULL) free(gradrow);
607	+
608	+	if (ra != NULL)                 /* extract eigenvectors & radii */
609	+	eigenvectors(uv, ra, hessian);
610	+	if (pg != NULL) {               /* tangential position gradient */
611	+	pg[0] = DOT(gradient, uv[0]);
612	+	pg[1] = DOT(gradient, uv[1]);
613	+	}
614	+	}
615	+
616	+
617	+	/* Compute direction gradient from a hemispherical sampling */
618	+	static void
619	+	ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
620	+	{
621	+	AMBSAMP *ap;
622	+	double  dgsum[2];
623	+	int     n;
624	+	FVECT   vd;
625	+	double  gfact;
626	+
627	+	dgsum[0] = dgsum[1] = 0.0;      /* sum values times -tan(theta) */
628	+	for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
629	+	/* use vector for azimuth + 90deg */
630	+	VSUB(vd, ap->p, hp->rp->rop);
631	+	/* brightness over cosine factor */
632	+	gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
633	+	/* sine = proj_radius/vd_length */
634	+	dgsum[0] -= DOT(uv[1], vd) * gfact;
635	+	dgsum[1] += DOT(uv[0], vd) * gfact;
636	+	}
637	+	dg[0] = dgsum[0] / (hp->ns*hp->ns);
638	+	dg[1] = dgsum[1] / (hp->ns*hp->ns);
639	+	}
640	+
641	+
642	+	/* Compute potential light leak direction flags for cache value */
643	+	static uint32
644	+	ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
645	+	{
646	+	const double    max_d = 1.0/(minarad*ambacc + 0.001);
647	+	const double    ang_res = 0.5*PI/hp->ns;
648	+	const double    ang_step = ang_res/((int)(16/PI*ang_res) + 1.01);
649	+	double          avg_d = 0;
650	+	uint32          flgs = 0;
651	+	FVECT           vec;
652	+	double          u, v;
653	+	double          ang, a1;
654	+	int             i, j;
655	+	/* don't bother for a few samples */
656	+	if (hp->ns < 8)
657	+	return(0);
658	+	/* check distances overhead */
659	+	for (i = hp->ns*3/4; i-- > hp->ns>>2; )
660	+	for (j = hp->ns*3/4; j-- > hp->ns>>2; )
661	+	avg_d += ambsam(hp,i,j).d;
662	+	avg_d *= 4.0/(hp->ns*hp->ns);
663	+	if (avg_d*r0 >= 1.0)            /* ceiling too low for corral? */
664	+	return(0);
665	+	if (avg_d >= max_d)             /* insurance */
666	+	return(0);
667	+	/* else circle around perimeter */
668	+	for (i = 0; i < hp->ns; i++)
669	+	for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) {
670	+	AMBSAMP *ap = &ambsam(hp,i,j);
671	+	if ((ap->d <= FTINY) | (ap->d >= max_d))
672	+	continue;       /* too far or too near */
673	+	VSUB(vec, ap->p, hp->rp->rop);
674	+	u = DOT(vec, uv[0]);
675	+	v = DOT(vec, uv[1]);
676	+	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v)
677	+	continue;       /* occluder outside ellipse */
678	+	ang = atan2a(v, u);     /* else set direction flags */
679	+	for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step)
680	+	flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0)));
681	+	}
682	+	return(flgs);
683	+	}
684	+
685	+
686	+	int
687	+	doambient(                              /* compute ambient component */
688	+	COLOR   rcol,                   /* input/output color */
689	+	RAY     *r,
690	+	double  wt,
691	+	FVECT   uv[2],                  /* returned (optional) */
692	+	float   ra[2],                  /* returned (optional) */
693	+	float   pg[2],                  /* returned (optional) */
694	+	float   dg[2],                  /* returned (optional) */
695	+	uint32  *crlp                   /* returned (optional) */
696	+	)
697	+	{
698	+	AMBHEMI *hp = samp_hemi(rcol, r, wt);
699	+	FVECT   my_uv[2];
700	+	double  d, K;
701	+	AMBSAMP *ap;
702	+	int     i;
703	+	/* clear return values */
704	+	if (uv != NULL)
705	+	memset(uv, 0, sizeof(FVECT)*2);
706	+	if (ra != NULL)
707	+	ra[0] = ra[1] = 0.0;
708	+	if (pg != NULL)
709	+	pg[0] = pg[1] = 0.0;
710	+	if (dg != NULL)
711	+	dg[0] = dg[1] = 0.0;
712	+	if (crlp != NULL)
713	+	*crlp = 0;
714	+	if (hp == NULL)                 /* sampling falure? */
715	+	return(0);
716	+
717	+	if ((ra == NULL) & (pg == NULL) & (dg == NULL) ||
718	+	(hp->sampOK < 0) | (hp->ns < 6)) {
719	+	free(hp);               /* Hessian not requested/possible */
720	+	return(-1);             /* value-only return value */
721	+	}
722	+	if ((d = bright(rcol)) > FTINY) {       /* normalize Y values */
723	+	d = 0.99*(hp->ns*hp->ns)/d;
724	+	K = 0.01;
725	+	} else {                        /* or fall back on geometric Hessian */
726	+	K = 1.0;
727	+	pg = NULL;
728	+	dg = NULL;
729	+	crlp = NULL;
730	+	}
731	+	ap = hp->sa;                    /* relative Y channel from here on... */
732	+	for (i = hp->ns*hp->ns; i--; ap++)
733	+	colval(ap->v,CIEY) = bright(ap->v)*d + K;
734	+
735	+	if (uv == NULL)                 /* make sure we have axis pointers */
736	+	uv = my_uv;
737	+	/* compute radii & pos. gradient */
738	+	ambHessian(hp, uv, ra, pg);
739	+
740	+	if (dg != NULL)                 /* compute direction gradient */
741	+	ambdirgrad(hp, uv, dg);
742	+
743	+	if (ra != NULL) {               /* scale/clamp radii */
744	+	if (pg != NULL) {
745	+	if (ra[0]*(d = fabs(pg[0])) > 1.0)
746	+	ra[0] = 1.0/d;
747	+	if (ra[1]*(d = fabs(pg[1])) > 1.0)
748	+	ra[1] = 1.0/d;
749	+	if (ra[0] > ra[1])
750	+	ra[0] = ra[1];
751	+	}
752	+	if (ra[0] < minarad) {
753	+	ra[0] = minarad;
754	+	if (ra[1] < minarad)
755	+	ra[1] = minarad;
756	+	}
757	+	ra[0] *= d = 1.0/sqrt(wt);
758	+	if ((ra[1] *= d) > 2.0*ra[0])
759	+	ra[1] = 2.0*ra[0];
760	+	if (ra[1] > maxarad) {
761	+	ra[1] = maxarad;
762	+	if (ra[0] > maxarad)
763	+	ra[0] = maxarad;
764	+	}
765	+	/* flag encroached directions */
766	+	if (crlp != NULL)
767	+	*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc);
768	+	if (pg != NULL) {       /* cap gradient if necessary */
769	+	d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];
770	+	if (d > 1.0) {
771	+	d = 1.0/sqrt(d);
772	+	pg[0] *= d;
773	+	pg[1] *= d;
774	+	}
775	+	}
776	+	}
777	+	free(hp);                       /* clean up and return */
778	+	return(1);
779	+	}
780	+
781	+
782	+	#else /* ! NEWAMB */
783	+
784	+
785		void
786		inithemi(                       /* initialize sampling hemisphere */
787	<	register AMBHEMI  *hp,
787	>	AMBHEMI  *hp,
788		COLOR ac,
789		RAY  *r,
790		double  wt
791		)
792		{
793		double  d;
794	<	register int  i;
794	>	int  i;
795		/* set number of divisions */
796		if (ambacc <= FTINY &&
797		wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
#	Line 58 | Line 824 | inithemi(                      /* initialize sampling hemisphere */
824
825		int
826		divsample(                              /* sample a division */
827	<	register AMBSAMP  *dp,
827	>	AMBSAMP  *dp,
828		AMBHEMI  *h,
829		RAY  *r
830		)
#	Line 69 | Line 835 | divsample(                             /* sample a division */
835		double  xd, yd, zd;
836		double  b2;
837		double  phi;
838	<	register int  i;
838	>	int  i;
839		/* ambient coefficient for weight */
840		if (ambacc > FTINY)
841		setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
#	Line 94 | Line 860 | divsample(                             /* sample a division */
860		ar.rdir[i] =    xd*h->ux[i] +
861		yd*h->uy[i] +
862		zd*h->uz[i];
863	+	checknorm(ar.rdir);
864		dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
865		rayvalue(&ar);
866		ndims--;
867		multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
868		addcolor(dp->v, ar.rcol);
869	<	/* use rt to improve gradient calc */
870	<	if (ar.rt > FTINY && ar.rt < FHUGE)
871	<	dp->r += 1.0/ar.rt;
869	>	/* use rxt to improve gradient calc */
870	>	if (ar.rxt > FTINY && ar.rxt < FHUGE)
871	>	dp->r += 1.0/ar.rxt;
872		/* (re)initialize error */
873		if (dp->n++) {
874		b2 = bright(dp->v)/dp->n - bright(ar.rcol);
#	Line 138 | Line 905 | ambnorm(                               /* standard order */
905		{
906		const AMBSAMP   *d1 = (const AMBSAMP *)p1;
907		const AMBSAMP   *d2 = (const AMBSAMP *)p2;
908	<	register int    c;
908	>	int     c;
909
910		if ( (c = d1->t - d2->t) )
911		return(c);
#	Line 148 | Line 915 | ambnorm(                               /* standard order */
915
916		double
917		doambient(                              /* compute ambient component */
918	<	COLOR  acol,
918	>	COLOR  rcol,
919		RAY  *r,
920		double  wt,
921		FVECT  pg,
922		FVECT  dg
923		)
924		{
925	<	double  b, d;
925	>	double  b, d=0;
926		AMBHEMI  hemi;
927		AMBSAMP  *div;
928		AMBSAMP  dnew;
929	<	register AMBSAMP  *dp;
929	>	double  acol[3];
930	>	AMBSAMP  *dp;
931		double  arad;
932		int  divcnt;
933	<	register int  i, j;
933	>	int  i, j;
934		/* initialize hemisphere */
935	<	inithemi(&hemi, acol, r, wt);
935	>	inithemi(&hemi, rcol, r, wt);
936		divcnt = hemi.nt * hemi.np;
937		/* initialize */
938		if (pg != NULL)
939		pg[0] = pg[1] = pg[2] = 0.0;
940		if (dg != NULL)
941		dg[0] = dg[1] = dg[2] = 0.0;
942	<	setcolor(acol, 0.0, 0.0, 0.0);
942	>	setcolor(rcol, 0.0, 0.0, 0.0);
943		if (divcnt == 0)
944		return(0.0);
945		/* allocate super-samples */
#	Line 183 | Line 951 | doambient(                             /* compute ambient component */
951		div = NULL;
952		/* sample the divisions */
953		arad = 0.0;
954	+	acol[0] = acol[1] = acol[2] = 0.0;
955		if ((dp = div) == NULL)
956		dp = &dnew;
957		divcnt = 0;
#	Line 204 | Line 973 | doambient(                             /* compute ambient component */
973		else
974		addcolor(acol, dp->v);
975		}
976	<	if (!divcnt)
976	>	if (!divcnt) {
977	>	if (div != NULL)
978	>	free((void *)div);
979		return(0.0);            /* no samples taken */
980	+	}
981		if (divcnt < hemi.nt*hemi.np) {
982		pg = dg = NULL;         /* incomplete sampling */
983		hemi.ns = 0;
#	Line 261 | Line 1033 | doambient(                             /* compute ambient component */
1033		}
1034		free((void *)div);
1035		}
1036	+	copycolor(rcol, acol);
1037		if (arad <= FTINY)
1038		arad = maxarad;
1039		else
#	Line 287 | Line 1060 | doambient(                             /* compute ambient component */
1060		void
1061		comperrs(                       /* compute initial error estimates */
1062		AMBSAMP  *da,   /* assumes standard ordering */
1063	<	register AMBHEMI  *hp
1063	>	AMBHEMI  *hp
1064		)
1065		{
1066		double  b, b2;
1067		int  i, j;
1068	<	register AMBSAMP  *dp;
1068	>	AMBSAMP  *dp;
1069		/* sum differences from neighbors */
1070		dp = da;
1071		for (i = 0; i < hp->nt; i++)
#	Line 340 | Line 1113 | void
1113		posgradient(                                    /* compute position gradient */
1114		FVECT  gv,
1115		AMBSAMP  *da,                   /* assumes standard ordering */
1116	<	register AMBHEMI  *hp
1116	>	AMBHEMI  *hp
1117		)
1118		{
1119	<	register int  i, j;
1119	>	int  i, j;
1120		double  nextsine, lastsine, b, d;
1121		double  mag0, mag1;
1122		double  phi, cosp, sinp, xd, yd;
1123	<	register AMBSAMP  *dp;
1123	>	AMBSAMP  *dp;
1124
1125		xd = yd = 0.0;
1126		for (j = 0; j < hp->np; j++) {
#	Line 398 | Line 1171 | void
1171		dirgradient(                                    /* compute direction gradient */
1172		FVECT  gv,
1173		AMBSAMP  *da,                   /* assumes standard ordering */
1174	<	register AMBHEMI  *hp
1174	>	AMBHEMI  *hp
1175		)
1176		{
1177	<	register int  i, j;
1177	>	int  i, j;
1178		double  mag;
1179		double  phi, xd, yd;
1180	<	register AMBSAMP  *dp;
1180	>	AMBSAMP  *dp;
1181
1182		xd = yd = 0.0;
1183		for (j = 0; j < hp->np; j++) {
#	Line 427 | Line 1200 | dirgradient(                                   /* compute direction gradient */
1200		for (i = 0; i < 3; i++)
1201		gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
1202		}
1203	+
1204	+	#endif  /* ! NEWAMB */

Diff Legend

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