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

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.55 by greg, Fri May 9 16:05:09 2014 UTC vs.
Revision 2.90 by greg, Wed Nov 15 18:02:52 2023 UTC

#	Line 21 \| Line 21 \| static const char RCSid[] = "$Id$";
21		#include "ambient.h"
22		#include "random.h"
23
24	<	#ifdef NEWAMB
24	>	#ifndef MINADIV
25	>	#define MINADIV 7 /* minimum # divisions in each dimension */
26	>	#endif
27
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	–
28		typedef struct {
52	–	COLOR v; /* hemisphere sample value */
53	–	float d; /* reciprocal distance (1/rt) */
29		FVECT p; /* intersection point */
30	+	float d; /* reciprocal distance */
31	+	SCOLOR v; /* hemisphere sample value */
32		} AMBSAMP; /* sample value */
33
34		typedef struct {
35		RAY rp; / originating ray sample */
59	–	FVECT ux, uy; /* tangent axis unit vectors */
36		int ns; /* number of samples per axis */
37	<	COLOR acoef; /* division contribution coefficient */
37	>	int sampOK; /* acquired full sample set? */
38	>	SCOLOR acoef; /* division contribution coefficient */
39	>	SCOLOR acol; /* 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 ambndx(h,i,j) ((i)*(h)->ns + (j))
45	<	#define ambsam(h,i,j) (h)->sa[ambndx(h,i,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, rcp, rI2_eJ2;
49		double I1, I2;
71	–	int valid;
50		} FFTRI; /* vectors and coefficients for Hessian calculation */
51
52
75	–	/* Get index for adjacent vertex */
53		static int
54	<	adjacent_verti(AMBHEMI *hp, int i, int j, int dbit)
54	>	ambcollision( /* proposed direciton collides? */
55	>	AMBHEMI *hp,
56	>	int i,
57	>	int j,
58	>	FVECT dv
59	>	)
60		{
61	<	int i0 = i*hp->ns + j;
62	<
63	<	switch (dbit) {
64	<	case VDB_y: return(i0 - hp->ns);
65	<	case VDB_x: return(i0 - 1);
66	<	case VDB_Xy: return(i0 - hp->ns + 1);
67	<	case VDB_xY: return(i0 + hp->ns - 1);
68	<	case VDB_X: return(i0 + 1);
69	<	case VDB_Y: return(i0 + hp->ns);
70	<	/* the following should never occur */
71	<	case VDB_xy: return(i0 - hp->ns - 1);
72	<	case VDB_XY: return(i0 + hp->ns + 1);
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. - .5cos_threshcos_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(-1);
86	>	return(0); /* nothing to worry about */
87		}
88
89
96	–	/* Get vertex direction bit for the opposite edge to complete triangle */
90		static int
91	<	vdb_edge(int db1, int db2)
92	<	{
93	<	switch (db1) {
94	<	case VDB_x: return(db2==VDB_y ? VDB_Xy : VDB_Y);
95	<	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(CONSISTENCY, "forbidden diagonal in vdb_edge()");
109	<	return(-1);
110	<	}
111	<
112	<
113	<	static AMBHEMI *
114	<	inithemi( /* initialize sampling hemisphere */
115	<	COLOR ac,
116	<	RAY *r,
117	<	double wt
91	>	ambsample( /* initial ambient division sample */
92	>	AMBHEMI *hp,
93	>	int i,
94	>	int j,
95	>	int n
96		)
97		{
98	<	AMBHEMI *hp;
99	<	double d;
122	<	int n, i;
123	<	/* set number of divisions */
124	<	if (ambacc <= FTINY &&
125	<	wt > (d = 0.8intens(ac)r->rweight/(ambdiv*minweight)))
126	<	wt = d; /* avoid ray termination */
127	<	n = sqrt(ambdiv * wt) + 0.5;
128	<	i = 1 + 5(ambacc > FTINY); / minimum number of samples */
129	<	if (n < i)
130	<	n = i;
131	<	/* allocate sampling array */
132	<	hp = (AMBHEMI )malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)(n*n - 1));
133	<	if (hp == NULL)
134	<	return(NULL);
135	<	hp->rp = r;
136	<	hp->ns = n;
137	<	/* assign coefficient */
138	<	copycolor(hp->acoef, ac);
139	<	d = 1.0/(n*n);
140	<	scalecolor(hp->acoef, d);
141	<	/* make tangent plane axes */
142	<	hp->uy[0] = 0.5 - frandom();
143	<	hp->uy[1] = 0.5 - frandom();
144	<	hp->uy[2] = 0.5 - frandom();
145	<	for (i = 3; i--; )
146	<	if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
147	<	break;
148	<	if (i < 0)
149	<	error(CONSISTENCY, "bad ray direction in inithemi");
150	<	hp->uy[i] = 1.0;
151	<	VCROSS(hp->ux, hp->uy, r->ron);
152	<	normalize(hp->ux);
153	<	VCROSS(hp->uy, r->ron, hp->ux);
154	<	/* we're ready to sample */
155	<	return(hp);
156	<	}
157	<
158	<
159	<	/* Sample ambient division and apply weighting coefficient */
160	<	static int
161	<	getambsamp(RAY arp, AMBHEMI hp, int i, int j, int n)
162	<	{
98	>	AMBSAMP *ap = &ambsam(hp,i,j);
99	>	RAY ar;
100		int hlist[3], ii;
101	<	double spt[2], zd;
101	>	RREAL spt[2];
102	>	double zd;
103	>	/* generate hemispherical sample */
104		/* ambient coefficient for weight */
105		if (ambacc > FTINY)
106	<	setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL);
106	>	setscolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
107		else
108	<	copycolor(arp->rcoef, hp->acoef);
109	<	if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0)
108	>	copyscolor(ar.rcoef, hp->acoef);
109	>	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
110		return(0);
111		if (ambacc > FTINY) {
112	<	multcolor(arp->rcoef, hp->acoef);
113	<	scalecolor(arp->rcoef, 1./AVGREFL);
112	>	smultscolor(ar.rcoef, hp->acoef);
113	>	scalescolor(ar.rcoef, 1./AVGREFL);
114		}
115		hlist[0] = hp->rp->rno;
116		hlist[1] = j;
117		hlist[2] = i;
118		multisamp(spt, 2, urand(ilhash(hlist,3)+n));
119	<	if (!n) { /* avoid border samples for n==0 */
120	<	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))
184	<	spt[1] = 0.1 + 0.8*frandom();
185	<	}
186	<	SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
119	>	resample:
120	>	square2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
121		zd = sqrt(1. - spt[0]spt[0] - spt[1]spt[1]);
122		for (ii = 3; ii--; )
123	<	arp->rdir[ii] = spt[0]*hp->ux[ii] +
123	>	ar.rdir[ii] = spt[0]*hp->ux[ii] +
124		spt[1]*hp->uy[ii] +
125		zd*hp->rp->ron[ii];
126	<	checknorm(arp->rdir);
127	<	dimlist[ndims++] = ambndx(hp,i,j) + 90171;
128	<	rayvalue(arp); /* evaluate ray */
129	<	ndims--; /* apply coefficient */
130	<	multcolor(arp->rcol, arp->rcoef);
126	>	checknorm(ar.rdir);
127	>	/* avoid coincident samples */
128	>	if (!n && ambcollision(hp, i, j, ar.rdir)) {
129	>	spt[0] = frandom(); spt[1] = frandom();
130	>	goto resample; /* reject this sample */
131	>	}
132	>	dimlist[ndims++] = AI(hp,i,j) + 90171;
133	>	rayvalue(&ar); /* evaluate ray */
134	>	ndims--;
135	>	zd = raydistance(&ar);
136	>	if (zd <= FTINY)
137	>	return(0); /* should never happen */
138	>	smultscolor(ar.rcol, ar.rcoef); /* apply coefficient */
139	>	if (zdap->d < 1.0) / new/closer distance? */
140	>	ap->d = 1.0/zd;
141	>	if (!n) { /* record first vertex & value */
142	>	if (zd > 10.0*thescene.cusize + 1000.)
143	>	zd = 10.0*thescene.cusize + 1000.;
144	>	VSUM(ap->p, ar.rorg, ar.rdir, zd);
145	>	copyscolor(ap->v, ar.rcol);
146	>	} else { /* else update recorded value */
147	>	sopscolor(hp->acol, -=, ap->v);
148	>	zd = 1.0/(double)(n+1);
149	>	scalescolor(ar.rcol, zd);
150	>	zd *= (double)n;
151	>	scalescolor(ap->v, zd);
152	>	saddscolor(ap->v, ar.rcol);
153	>	}
154	>	saddscolor(hp->acol, ap->v); /* add to our sum */
155		return(1);
156		}
157
158
159	<	static AMBSAMP *
202	<	ambsample( /* initial ambient division sample */
203	<	AMBHEMI *hp,
204	<	int i,
205	<	int j
206	<	)
207	<	{
208	<	AMBSAMP *ap = &ambsam(hp,i,j);
209	<	RAY ar;
210	<	/* generate hemispherical sample */
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;
218	<	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
219	<	copycolor(ap->v, ar.rcol);
220	<	return(ap);
221	<	}
222	<
223	<
224	<	/* Estimate errors based on ambient division differences */
159	>	/* Estimate variance based on ambient division differences */
160		static float *
161		getambdiffs(AMBHEMI *hp)
162		{
163	+	const double normf = 1./bright(hp->acoef);
164		float earr = (float )calloc(hp->ns*hp->ns, sizeof(float));
165		float *ep;
166		AMBSAMP *ap;
167	<	double b, d2;
167	>	double b, b1, d2;
168		int i, j;
169
170		if (earr == NULL) /* out of memory? */
171		return(NULL);
172	<	/* compute squared neighbor diffs */
172	>	/* sum squared neighbor diffs */
173		for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
174		for (j = 0; j < hp->ns; j++, ap++, ep++) {
175	<	b = bright(ap[0].v);
175	>	b = pbright(ap[0].v);
176		if (i) { /* from above */
177	<	d2 = b - bright(ap[-hp->ns].v);
178	<	d2 *= d2;
177	>	b1 = pbright(ap[-hp->ns].v);
178	>	d2 = b - b1;
179	>	d2 = d2normf/(b + b1 + FTINY);
180		ep[0] += d2;
181		ep[-hp->ns] += d2;
182		}
183		if (!j) continue;
184		/* from behind */
185	<	d2 = b - bright(ap[-1].v);
186	<	d2 *= d2;
185	>	b1 = pbright(ap[-1].v);
186	>	d2 = b - b1;
187	>	d2 = d2normf/(b + b1 + FTINY);
188		ep[0] += d2;
189		ep[-1] += d2;
190		if (!i) continue;
191		/* diagonal */
192	<	d2 = b - bright(ap[-hp->ns-1].v);
193	<	d2 *= d2;
192	>	b1 = pbright(ap[-hp->ns-1].v);
193	>	d2 = b - b1;
194	>	d2 = d2normf/(b + b1 + FTINY);
195		ep[0] += d2;
196		ep[-hp->ns-1] += d2;
197		}
#	Line 275 \| Line 214 \| *getambdiffs(AMBHEMI hp)**
214
215		/* Perform super-sampling on hemisphere (introduces bias) */
216		static void
217	<	ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
217	>	ambsupersamp(AMBHEMI *hp, int cnt)
218		{
219		float *earr = getambdiffs(hp);
220		double e2rem = 0;
282	–	AMBSAMP *ap;
283	–	RAY ar;
284	–	double asum[3];
221		float *ep;
222		int i, j, n, nss;
223
#	Line 291 \| Line 227 \| *ambsupersamp(double acol[3], AMBHEMI hp, int cnt)**
227		for (ep = earr + hp->ns*hp->ns; ep > earr; )
228		e2rem += *--ep;
229		ep = earr; /* perform super-sampling */
230	<	for (ap = hp->sa, i = 0; i < hp->ns; i++)
231	<	for (j = 0; j < hp->ns; j++, ap++) {
230	>	for (i = 0; i < hp->ns; i++)
231	>	for (j = 0; j < hp->ns; j++) {
232		if (e2rem <= FTINY)
233		goto done; /* nothing left to do */
234		nss = ep/e2remcnt + frandom();
235	<	asum[0] = asum[1] = asum[2] = 0.0;
236	<	for (n = 1; n <= nss; n++) {
237	<	if (!getambsamp(&ar, hp, i, j, n)) {
302	<	nss = n-1;
303	<	break;
304	<	}
305	<	addcolor(asum, ar.rcol);
306	<	}
307	<	if (nss) { /* update returned ambient value */
308	<	const double ssf = 1./(nss + 1.);
309	<	for (n = 3; n--; )
310	<	acol[n] += ssf*asum[n] +
311	<	(ssf - 1.)*colval(ap->v,n);
312	<	}
313	<	e2rem -= ep++; / update remainders */
314	<	cnt -= nss;
235	>	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
236	>	if (!--cnt) goto done;
237	>	e2rem -= ep++; / update remainder */
238		}
239		done:
240		free(earr);
241		}
242
243
244	<	/* Compute vertex flags, indicating farthest in each direction */
245	<	static uby8 *
246	<	vertex_flags(AMBHEMI *hp)
244	>	static AMBHEMI *
245	>	samp_hemi( /* sample indirect hemisphere */
246	>	SCOLOR rcol,
247	>	RAY *r,
248	>	double wt
249	>	)
250		{
251	<	uby8 vflags = (uby8 )calloc(hp->ns*hp->ns, sizeof(uby8));
252	<	uby8 *vf;
253	<	AMBSAMP *ap;
254	<	int i, j;
255	<
256	<	if (vflags == NULL)
257	<	error(SYSTEM, "out of memory in vertex_flags()");
258	<	vf = vflags;
259	<	ap = hp->sa; /* compute farthest along first row */
260	<	for (j = 0; j < hp->ns-1; j++, vf++, ap++)
261	<	if (ap[0].d <= ap[1].d)
262	<	vf[0] \|= 1<<VDB_X;
263	<	else
264	<	vf[1] \|= 1<<VDB_x;
265	<	++vf; ++ap;
266	<	/* flag subsequent rows */
267	<	for (i = 1; i < hp->ns; i++) {
268	<	for (j = 0; j < hp->ns-1; j++, vf++, ap++) {
269	<	if (ap[0].d <= ap[-hp->ns].d) /* row before */
270	<	vf[0] \|= 1<<VDB_y;
271	<	else
272	<	vf[-hp->ns] \|= 1<<VDB_Y;
273	<	if (ap[0].d <= ap[1-hp->ns].d) /* diagonal we care about */
274	<	vf[0] \|= 1<<VDB_Xy;
275	<	else
276	<	vf[1-hp->ns] \|= 1<<VDB_xY;
277	<	if (ap[0].d <= ap[1].d) /* column after */
278	<	vf[0] \|= 1<<VDB_X;
279	<	else
280	<	vf[1] \|= 1<<VDB_x;
281	<	}
282	<	if (ap[0].d <= ap[-hp->ns].d) /* final column edge */
283	<	vf[0] \|= 1<<VDB_y;
284	<	else
285	<	vf[-hp->ns] \|= 1<<VDB_Y;
286	<	++vf; ++ap;
251	>	AMBHEMI *hp;
252	>	double d;
253	>	int n, i, j;
254	>	/* insignificance check */
255	>	d = sintens(rcol);
256	>	if (d <= FTINY)
257	>	return(NULL);
258	>	/* set number of divisions */
259	>	if (ambacc <= FTINY &&
260	>	wt > (d = 0.8r->rweight/(ambdiv*minweight)))
261	>	wt = d; /* avoid ray termination */
262	>	n = sqrt(ambdiv * wt) + 0.5;
263	>	i = 1 + (MINADIV-1)*(ambacc > FTINY);
264	>	if (n < i) /* use minimum number of samples? */
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	>	scolorblack(hp->acol);
273	>	memset(hp->sa, 0, sizeof(AMBSAMP)nn);
274	>	hp->sampOK = 0;
275	>	/* assign coefficient */
276	>	copyscolor(hp->acoef, rcol);
277	>	d = 1.0/(n*n);
278	>	scalescolor(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	>	copyscolor(rcol, hp->acol);
288	>	if (!hp->sampOK) { /* utter failure? */
289	>	free(hp);
290	>	return(NULL);
291		}
292	<	return(vflags);
292	>	if (hp->sampOK < hp->ns*hp->ns) {
293	>	hp->sampOK = -1; / soft failure */
294	>	return(hp);
295	>	}
296	>	if (hp->sampOK <= MINADIV*MINADIV)
297	>	return(hp); /* don't bother super-sampling */
298	>	n = ambssamp*wt + 0.5;
299	>	if (n > 8) { /* perform super-sampling? */
300	>	ambsupersamp(hp, n);
301	>	copyscolor(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, int i, int j, int dbit1, int dbit2, const uby8 vflags)
309	>	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
310		{
311	<	const int v0 = ambndx(hp,i,j);
312	<	const int tflags = (1<<dbit1 \| 1<<dbit2);
313	<	int v1, v2;
314	<
315	<	if ((vflags[v0] & tflags) == tflags) /* is v0 the farthest? */
316	<	return(colval(hp->sa[v0].v,CIEY));
317	<	v1 = adjacent_verti(hp, i, j, dbit1);
318	<	if (vflags[v0] & 1<<dbit2) /* v1 farthest if v0>v2 */
378	<	return(colval(hp->sa[v1].v,CIEY));
379	<	v2 = adjacent_verti(hp, i, j, dbit2);
380	<	if (vflags[v0] & 1<<dbit1) /* v2 farthest if v0>v1 */
381	<	return(colval(hp->sa[v2].v,CIEY));
382	<	/* else check if v1>v2 */
383	<	if (vflags[v1] & 1<<vdb_edge(dbit1,dbit2))
384	<	return(colval(hp->sa[v1].v,CIEY));
385	<	return(colval(hp->sa[v2].v,CIEY));
311	>	if (hp->sa[n1].d <= hp->sa[n2].d) {
312	>	if (hp->sa[n1].d <= hp->sa[n3].d)
313	>	return(hp->sa[n1].v[0]);
314	>	return(hp->sa[n3].v[0]);
315	>	}
316	>	if (hp->sa[n2].d <= hp->sa[n3].d)
317	>	return(hp->sa[n2].v[0]);
318	>	return(hp->sa[n3].v[0]);
319		}
320
321
322		/* Compute vectors and coefficients for Hessian/gradient calcs */
323		static void
324	<	comp_fftri(FFTRI ftp, AMBHEMI hp, int i, int j, int dbit, const uby8 *vflags)
324	>	comp_fftri(FFTRI ftp, AMBHEMI hp, const int n0, const int n1)
325		{
326	<	const int i0 = ambndx(hp,i,j);
327	<	double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
395	<	int i1, ii;
326	>	double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
327	>	int ii;
328
329	<	ftp->valid = 0; /* check if we can skip this edge */
330	<	ii = adjacent_trifl[dbit];
331	<	if ((vflags[i0] & ii) == ii) /* cancels if vertex used as value */
400	<	return;
401	<	i1 = adjacent_verti(hp, i, j, dbit);
402	<	ii = adjacent_trifl[VDB_OPP(dbit)];
403	<	if ((vflags[i1] & ii) == ii) /* on either end (for both triangles) */
404	<	return;
405	<	/* else go ahead with calculation */
406	<	VSUB(ftp->r_i, hp->sa[i0].p, hp->rp->rop);
407	<	VSUB(ftp->r_i1, hp->sa[i1].p, hp->rp->rop);
408	<	VSUB(ftp->e_i, hp->sa[i1].p, hp->sa[i0].p);
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);
#	Line 419 \| Line 342 \| *comp_fftri(FFTRI ftp, AMBHEMI hp, int i, int j, int*
342		J2 = ( 0.5(rdot_r - rdot_r1) - dot_erftp->I2 ) / dot_e;
343		for (ii = 3; ii--; )
344		ftp->rI2_eJ2[ii] = ftp->I2ftp->r_i[ii] + J2ftp->e_i[ii];
422	–	ftp->valid++;
345		}
346
347
#	Line 445 \| Line 367 \| *comp_hessian(FVECT hess[3], FFTRI ftp, FVECT nrm)**
367		double d1, d2, d3, d4;
368		double I3, J3, K3;
369		int i, j;
448	–
449	–	if (!ftp->valid) { /* preemptive test */
450	–	memset(hess, 0, sizeof(FVECT)*3);
451	–	return;
452	–	}
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);
#	Line 513 \| Line 430 \| *comp_gradient(FVECT grad, FFTRI ftp, FVECT nrm)**
430		double f1;
431		int i;
432
516	–	if (!ftp->valid) { /* preemptive test */
517	–	memset(grad, 0, sizeof(FVECT));
518	–	return;
519	–	}
433		f1 = 2.0*DOT(nrm, ftp->rcp);
434		VCROSS(ncp, nrm, ftp->e_i);
435		for (i = 3; i--; )
#	Line 606 \| Line 519 \| *ambHessian( / anisotropic radii & pos. gradient /*
519		static char memerrmsg[] = "out of memory in ambHessian()";
520		FVECT (*hessrow)[3] = NULL;
521		FVECT *gradrow = NULL;
609	–	uby8 *vflags;
522		FVECT hessian[3];
523		FVECT gradient;
524		FFTRI fftr;
#	Line 628 \| Line 540 \| *ambHessian( / anisotropic radii & pos. gradient /*
540		error(SYSTEM, memerrmsg);
541		memset(gradient, 0, sizeof(gradient));
542		}
631	–	/* get vertex position flags */
632	–	vflags = vertex_flags(hp);
543		/* compute first row of edges */
544		for (j = 0; j < hp->ns-1; j++) {
545	<	comp_fftri(&fftr, hp, 0, j, VDB_X, vflags);
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)
#	Line 642 \| Line 552 \| *ambHessian( / anisotropic radii & pos. gradient /*
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, i, 0, VDB_Y, vflags);
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)
#	Line 651 \| Line 561 \| *ambHessian( / anisotropic radii & pos. gradient /*
561		FVECT hessdia[3]; /* compute triangle contributions */
562		FVECT graddia;
563		double backg;
564	<	backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags);
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, i, j+1, VDB_xY, vflags);
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);
#	Line 665 \| Line 576 \| *ambHessian( / anisotropic radii & pos. gradient /*
576		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
577		}
578		/* initialize edge in next row */
579	<	comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags);
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, i+1, j+1, VDB_x, VDB_y, vflags);
586	<	comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags);
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);
#	Line 692 \| Line 604 \| *ambHessian( / anisotropic radii & pos. gradient /*
604		/* release row buffers */
605		if (hessrow != NULL) free(hessrow);
606		if (gradrow != NULL) free(gradrow);
695	–	free(vflags);
607
608		if (ra != NULL) /* extract eigenvectors & radii */
609		eigenvectors(uv, ra, hessian);
#	Line 718 \| Line 629 \| *ambdirgrad(AMBHEMI hp, FVECT uv[2], float dg[2])**
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);
632	>	gfact = ap->v[0] / 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;
#	Line 733 \| Line 644 \| 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-1);
648	<	const double ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY));
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 < 12)
656	>	if (hp->ns < 8)
657		return(0);
658		/* check distances overhead */
659		for (i = hp->ns*3/4; i-- > hp->ns>>2; )
#	Line 754 \| Line 668 \| *ambcorral(AMBHEMI hp, FVECT uv[2], const double r0, c**
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);
757	–	FVECT vec;
758	–	double u, v;
759	–	double ang, a1;
760	–	int abp;
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]) * ap->d;
675	<	v = DOT(vec, uv[1]) * ap->d;
676	<	if ((r0r0uu + r1r1vv) * ap->d*ap->d <= 1.0)
674	>	u = DOT(vec, uv[0]);
675	>	v = DOT(vec, uv[1]);
676	>	if ((r0r0uu + r1r1vv) * ap->dap->d <= uu + v*v)
677		continue; /* occluder outside ellipse */
678		ang = atan2a(v, u); /* else set direction flags */
679	<	for (a1 = ang-.5ang_res; a1 <= ang+.5ang_res; a1 += ang_step)
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);
#	Line 775 \| Line 685 \| *ambcorral(AMBHEMI hp, FVECT uv[2], const double r0, c**
685
686		int
687		doambient( /* compute ambient component */
688	<	COLOR rcol, /* input/output color */
688	>	SCOLOR rcol, /* input/output color */
689		RAY *r,
690		double wt,
691		FVECT uv[2], /* returned (optional) */
#	Line 785 \| Line 695 \| *doambient( / compute ambient component /*
695		uint32 crlp / returned (optional) */
696		)
697		{
698	<	AMBHEMI *hp = inithemi(rcol, r, wt);
789	<	int cnt;
698	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
699		FVECT my_uv[2];
700	<	double d, K, acol[3];
700	>	double d, K;
701		AMBSAMP *ap;
702	<	int i, j;
703	<	/* check/initialize */
795	<	if (hp == NULL)
796	<	return(0);
702	>	int i;
703	>	/* clear return values */
704		if (uv != NULL)
705		memset(uv, 0, sizeof(FVECT)*2);
706		if (ra != NULL)
#	Line 804 \| Line 711 \| *doambient( / compute ambient component /*
711		dg[0] = dg[1] = 0.0;
712		if (crlp != NULL)
713		*crlp = 0;
714	<	/* sample the hemisphere */
715	<	acol[0] = acol[1] = acol[2] = 0.0;
716	<	cnt = 0;
717	<	for (i = hp->ns; i--; )
718	<	for (j = hp->ns; j--; )
719	<	if ((ap = ambsample(hp, i, j)) != NULL) {
720	<	addcolor(acol, ap->v);
814	<	++cnt;
815	<	}
816	<	if (!cnt) {
817	<	setcolor(rcol, 0.0, 0.0, 0.0);
818	<	free(hp);
819	<	return(0); /* no valid samples */
714	>	if (hp == NULL) /* sampling falure? */
715	>	return(0);
716	>
717	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) \|\|
718	>	(hp->sampOK < 0) \| (hp->ns < MINADIV)) {
719	>	free(hp); /* Hessian not requested/possible */
720	>	return(-1); /* value-only return value */
721		}
722	<	if (cnt < hp->nshp->ns) { / incomplete sampling? */
723	<	copycolor(rcol, acol);
823	<	free(hp);
824	<	return(-1); /* return value w/o Hessian */
825	<	}
826	<	cnt = ambssampwt + 0.5; / perform super-sampling? */
827	<	if (cnt > 8)
828	<	ambsupersamp(acol, hp, cnt);
829	<	copycolor(rcol, acol); /* final indirect irradiance/PI */
830	<	if ((ra == NULL) & (pg == NULL) & (dg == NULL)) {
831	<	free(hp);
832	<	return(-1); /* no radius or gradient calc. */
833	<	}
834	<	if ((d = bright(acol)) > FTINY) { /* normalize Y values */
835	<	d = 0.99(hp->nshp->ns)/d;
722	>	if ((d = scolor_photopic(rcol)) > FTINY) {
723	>	d = 0.99(hp->nshp->ns)/d; /* normalize Y values */
724		K = 0.01;
725		} else { /* or fall back on geometric Hessian */
726		K = 1.0;
#	Line 840 \| Line 728 \| *doambient( / compute ambient component /*
728		dg = NULL;
729		crlp = NULL;
730		}
731	<	ap = hp->sa; /* relative Y channel from here on... */
731	>	ap = hp->sa; /* single channel from here on... */
732		for (i = hp->ns*hp->ns; i--; ap++)
733	<	colval(ap->v,CIEY) = bright(ap->v)*d + K;
733	>	ap->v[0] = scolor_mean(ap->v)*d + K;
734
735		if (uv == NULL) /* make sure we have axis pointers */
736		uv = my_uv;
#	Line 866 \| Line 754 \| *doambient( / compute ambient component /*
754		if (ra[1] < minarad)
755		ra[1] = minarad;
756		}
757	<	ra[0] *= d = 1.0/sqrt(sqrt(wt));
757	>	ra[0] *= d = 1.0/sqrt(wt);
758		if ((ra[1] = d) > 2.0ra[0])
759		ra[1] = 2.0*ra[0];
760		if (ra[1] > maxarad) {
#	Line 875 \| Line 763 \| *doambient( / compute ambient component /*
763		ra[0] = maxarad;
764		}
765		/* flag encroached directions */
766	<	if ((wt >= 0.5-FTINY) & (crlp != NULL))
766	>	if (crlp != NULL) /* XXX doesn't update with changes to ambacc */
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];
#	Line 889 \| Line 777 \| *doambient( / compute ambient component /*
777		free(hp); /* clean up and return */
778		return(1);
779		}
892	–
893	–
894	–	#else /* ! NEWAMB */
895	–
896	–
897	–	void
898	–	inithemi( /* initialize sampling hemisphere */
899	–	AMBHEMI *hp,
900	–	COLOR ac,
901	–	RAY *r,
902	–	double wt
903	–	)
904	–	{
905	–	double d;
906	–	int i;
907	–	/* set number of divisions */
908	–	if (ambacc <= FTINY &&
909	–	wt > (d = 0.8intens(ac)r->rweight/(ambdiv*minweight)))
910	–	wt = d; /* avoid ray termination */
911	–	hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
912	–	i = ambacc > FTINY ? 3 : 1; /* minimum number of samples */
913	–	if (hp->nt < i)
914	–	hp->nt = i;
915	–	hp->np = PI * hp->nt + 0.5;
916	–	/* set number of super-samples */
917	–	hp->ns = ambssamp * wt + 0.5;
918	–	/* assign coefficient */
919	–	copycolor(hp->acoef, ac);
920	–	d = 1.0/(hp->nt*hp->np);
921	–	scalecolor(hp->acoef, d);
922	–	/* make axes */
923	–	VCOPY(hp->uz, r->ron);
924	–	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
925	–	for (i = 0; i < 3; i++)
926	–	if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
927	–	break;
928	–	if (i >= 3)
929	–	error(CONSISTENCY, "bad ray direction in inithemi");
930	–	hp->uy[i] = 1.0;
931	–	fcross(hp->ux, hp->uy, hp->uz);
932	–	normalize(hp->ux);
933	–	fcross(hp->uy, hp->uz, hp->ux);
934	–	}
935	–
936	–
937	–	int
938	–	divsample( /* sample a division */
939	–	AMBSAMP *dp,
940	–	AMBHEMI *h,
941	–	RAY *r
942	–	)
943	–	{
944	–	RAY ar;
945	–	int hlist[3];
946	–	double spt[2];
947	–	double xd, yd, zd;
948	–	double b2;
949	–	double phi;
950	–	int i;
951	–	/* ambient coefficient for weight */
952	–	if (ambacc > FTINY)
953	–	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
954	–	else
955	–	copycolor(ar.rcoef, h->acoef);
956	–	if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
957	–	return(-1);
958	–	if (ambacc > FTINY) {
959	–	multcolor(ar.rcoef, h->acoef);
960	–	scalecolor(ar.rcoef, 1./AVGREFL);
961	–	}
962	–	hlist[0] = r->rno;
963	–	hlist[1] = dp->t;
964	–	hlist[2] = dp->p;
965	–	multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
966	–	zd = sqrt((dp->t + spt[0])/h->nt);
967	–	phi = 2.0PI (dp->p + spt[1])/h->np;
968	–	xd = tcos(phi) * zd;
969	–	yd = tsin(phi) * zd;
970	–	zd = sqrt(1.0 - zd*zd);
971	–	for (i = 0; i < 3; i++)
972	–	ar.rdir[i] = xd*h->ux[i] +
973	–	yd*h->uy[i] +
974	–	zd*h->uz[i];
975	–	checknorm(ar.rdir);
976	–	dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
977	–	rayvalue(&ar);
978	–	ndims--;
979	–	multcolor(ar.rcol, ar.rcoef); /* apply coefficient */
980	–	addcolor(dp->v, ar.rcol);
981	–	/* use rt to improve gradient calc */
982	–	if (ar.rt > FTINY && ar.rt < FHUGE)
983	–	dp->r += 1.0/ar.rt;
984	–	/* (re)initialize error */
985	–	if (dp->n++) {
986	–	b2 = bright(dp->v)/dp->n - bright(ar.rcol);
987	–	b2 = b2b2 + dp->k((dp->n-1)*(dp->n-1));
988	–	dp->k = b2/(dp->n*dp->n);
989	–	} else
990	–	dp->k = 0.0;
991	–	return(0);
992	–	}
993	–
994	–
995	–	static int
996	–	ambcmp( /* decreasing order */
997	–	const void *p1,
998	–	const void *p2
999	–	)
1000	–	{
1001	–	const AMBSAMP d1 = (const AMBSAMP )p1;
1002	–	const AMBSAMP d2 = (const AMBSAMP )p2;
1003	–
1004	–	if (d1->k < d2->k)
1005	–	return(1);
1006	–	if (d1->k > d2->k)
1007	–	return(-1);
1008	–	return(0);
1009	–	}
1010	–
1011	–
1012	–	static int
1013	–	ambnorm( /* standard order */
1014	–	const void *p1,
1015	–	const void *p2
1016	–	)
1017	–	{
1018	–	const AMBSAMP d1 = (const AMBSAMP )p1;
1019	–	const AMBSAMP d2 = (const AMBSAMP )p2;
1020	–	int c;
1021	–
1022	–	if ( (c = d1->t - d2->t) )
1023	–	return(c);
1024	–	return(d1->p - d2->p);
1025	–	}
1026	–
1027	–
1028	–	double
1029	–	doambient( /* compute ambient component */
1030	–	COLOR rcol,
1031	–	RAY *r,
1032	–	double wt,
1033	–	FVECT pg,
1034	–	FVECT dg
1035	–	)
1036	–	{
1037	–	double b, d=0;
1038	–	AMBHEMI hemi;
1039	–	AMBSAMP *div;
1040	–	AMBSAMP dnew;
1041	–	double acol[3];
1042	–	AMBSAMP *dp;
1043	–	double arad;
1044	–	int divcnt;
1045	–	int i, j;
1046	–	/* initialize hemisphere */
1047	–	inithemi(&hemi, rcol, r, wt);
1048	–	divcnt = hemi.nt * hemi.np;
1049	–	/* initialize */
1050	–	if (pg != NULL)
1051	–	pg[0] = pg[1] = pg[2] = 0.0;
1052	–	if (dg != NULL)
1053	–	dg[0] = dg[1] = dg[2] = 0.0;
1054	–	setcolor(rcol, 0.0, 0.0, 0.0);
1055	–	if (divcnt == 0)
1056	–	return(0.0);
1057	–	/* allocate super-samples */
1058	–	if (hemi.ns > 0 \|\| pg != NULL \|\| dg != NULL) {
1059	–	div = (AMBSAMP )malloc(divcntsizeof(AMBSAMP));
1060	–	if (div == NULL)
1061	–	error(SYSTEM, "out of memory in doambient");
1062	–	} else
1063	–	div = NULL;
1064	–	/* sample the divisions */
1065	–	arad = 0.0;
1066	–	acol[0] = acol[1] = acol[2] = 0.0;
1067	–	if ((dp = div) == NULL)
1068	–	dp = &dnew;
1069	–	divcnt = 0;
1070	–	for (i = 0; i < hemi.nt; i++)
1071	–	for (j = 0; j < hemi.np; j++) {
1072	–	dp->t = i; dp->p = j;
1073	–	setcolor(dp->v, 0.0, 0.0, 0.0);
1074	–	dp->r = 0.0;
1075	–	dp->n = 0;
1076	–	if (divsample(dp, &hemi, r) < 0) {
1077	–	if (div != NULL)
1078	–	dp++;
1079	–	continue;
1080	–	}
1081	–	arad += dp->r;
1082	–	divcnt++;
1083	–	if (div != NULL)
1084	–	dp++;
1085	–	else
1086	–	addcolor(acol, dp->v);
1087	–	}
1088	–	if (!divcnt) {
1089	–	if (div != NULL)
1090	–	free((void *)div);
1091	–	return(0.0); /* no samples taken */
1092	–	}
1093	–	if (divcnt < hemi.nt*hemi.np) {
1094	–	pg = dg = NULL; /* incomplete sampling */
1095	–	hemi.ns = 0;
1096	–	} else if (arad > FTINY && divcnt/arad < minarad) {
1097	–	hemi.ns = 0; /* close enough */
1098	–	} else if (hemi.ns > 0) { /* else perform super-sampling? */
1099	–	comperrs(div, &hemi); /* compute errors */
1100	–	qsort(div, divcnt, sizeof(AMBSAMP), ambcmp); /* sort divs */
1101	–	/* super-sample */
1102	–	for (i = hemi.ns; i > 0; i--) {
1103	–	dnew = *div;
1104	–	if (divsample(&dnew, &hemi, r) < 0) {
1105	–	dp++;
1106	–	continue;
1107	–	}
1108	–	dp = div; /* reinsert */
1109	–	j = divcnt < i ? divcnt : i;
1110	–	while (--j > 0 && dnew.k < dp[1].k) {
1111	–	dp = (dp+1);
1112	–	dp++;
1113	–	}
1114	–	*dp = dnew;
1115	–	}
1116	–	if (pg != NULL \|\| dg != NULL) /* restore order */
1117	–	qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
1118	–	}
1119	–	/* compute returned values */
1120	–	if (div != NULL) {
1121	–	arad = 0.0; /* note: divcnt may be < ntnp /
1122	–	for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
1123	–	arad += dp->r;
1124	–	if (dp->n > 1) {
1125	–	b = 1.0/dp->n;
1126	–	scalecolor(dp->v, b);
1127	–	dp->r *= b;
1128	–	dp->n = 1;
1129	–	}
1130	–	addcolor(acol, dp->v);
1131	–	}
1132	–	b = bright(acol);
1133	–	if (b > FTINY) {
1134	–	b = 1.0/b; /* compute & normalize gradient(s) */
1135	–	if (pg != NULL) {
1136	–	posgradient(pg, div, &hemi);
1137	–	for (i = 0; i < 3; i++)
1138	–	pg[i] *= b;
1139	–	}
1140	–	if (dg != NULL) {
1141	–	dirgradient(dg, div, &hemi);
1142	–	for (i = 0; i < 3; i++)
1143	–	dg[i] *= b;
1144	–	}
1145	–	}
1146	–	free((void *)div);
1147	–	}
1148	–	copycolor(rcol, acol);
1149	–	if (arad <= FTINY)
1150	–	arad = maxarad;
1151	–	else
1152	–	arad = (divcnt+hemi.ns)/arad;
1153	–	if (pg != NULL) { /* reduce radius if gradient large */
1154	–	d = DOT(pg,pg);
1155	–	if (daradarad > 1.0)
1156	–	arad = 1.0/sqrt(d);
1157	–	}
1158	–	if (arad < minarad) {
1159	–	arad = minarad;
1160	–	if (pg != NULL && daradarad > 1.0) { /* cap gradient */
1161	–	d = 1.0/arad/sqrt(d);
1162	–	for (i = 0; i < 3; i++)
1163	–	pg[i] *= d;
1164	–	}
1165	–	}
1166	–	if ((arad /= sqrt(wt)) > maxarad)
1167	–	arad = maxarad;
1168	–	return(arad);
1169	–	}
1170	–
1171	–
1172	–	void
1173	–	comperrs( /* compute initial error estimates */
1174	–	AMBSAMP da, / assumes standard ordering */
1175	–	AMBHEMI *hp
1176	–	)
1177	–	{
1178	–	double b, b2;
1179	–	int i, j;
1180	–	AMBSAMP *dp;
1181	–	/* sum differences from neighbors */
1182	–	dp = da;
1183	–	for (i = 0; i < hp->nt; i++)
1184	–	for (j = 0; j < hp->np; j++) {
1185	–	#ifdef DEBUG
1186	–	if (dp->t != i \|\| dp->p != j)
1187	–	error(CONSISTENCY,
1188	–	"division order in comperrs");
1189	–	#endif
1190	–	b = bright(dp[0].v);
1191	–	if (i > 0) { /* from above */
1192	–	b2 = bright(dp[-hp->np].v) - b;
1193	–	b2 = b2 0.25;
1194	–	dp[0].k += b2;
1195	–	dp[-hp->np].k += b2;
1196	–	}
1197	–	if (j > 0) { /* from behind */
1198	–	b2 = bright(dp[-1].v) - b;
1199	–	b2 = b2 0.25;
1200	–	dp[0].k += b2;
1201	–	dp[-1].k += b2;
1202	–	} else { /* around */
1203	–	b2 = bright(dp[hp->np-1].v) - b;
1204	–	b2 = b2 0.25;
1205	–	dp[0].k += b2;
1206	–	dp[hp->np-1].k += b2;
1207	–	}
1208	–	dp++;
1209	–	}
1210	–	/* divide by number of neighbors */
1211	–	dp = da;
1212	–	for (j = 0; j < hp->np; j++) /* top row */
1213	–	(dp++)->k *= 1.0/3.0;
1214	–	if (hp->nt < 2)
1215	–	return;
1216	–	for (i = 1; i < hp->nt-1; i++) /* central region */
1217	–	for (j = 0; j < hp->np; j++)
1218	–	(dp++)->k *= 0.25;
1219	–	for (j = 0; j < hp->np; j++) /* bottom row */
1220	–	(dp++)->k *= 1.0/3.0;
1221	–	}
1222	–
1223	–
1224	–	void
1225	–	posgradient( /* compute position gradient */
1226	–	FVECT gv,
1227	–	AMBSAMP da, / assumes standard ordering */
1228	–	AMBHEMI *hp
1229	–	)
1230	–	{
1231	–	int i, j;
1232	–	double nextsine, lastsine, b, d;
1233	–	double mag0, mag1;
1234	–	double phi, cosp, sinp, xd, yd;
1235	–	AMBSAMP *dp;
1236	–
1237	–	xd = yd = 0.0;
1238	–	for (j = 0; j < hp->np; j++) {
1239	–	dp = da + j;
1240	–	mag0 = mag1 = 0.0;
1241	–	lastsine = 0.0;
1242	–	for (i = 0; i < hp->nt; i++) {
1243	–	#ifdef DEBUG
1244	–	if (dp->t != i \|\| dp->p != j)
1245	–	error(CONSISTENCY,
1246	–	"division order in posgradient");
1247	–	#endif
1248	–	b = bright(dp->v);
1249	–	if (i > 0) {
1250	–	d = dp[-hp->np].r;
1251	–	if (dp[0].r > d) d = dp[0].r;
1252	–	/* sin(t)cos(t)^2 /
1253	–	d = lastsine (1.0 - (double)i/hp->nt);
1254	–	mag0 += d*(b - bright(dp[-hp->np].v));
1255	–	}
1256	–	nextsine = sqrt((double)(i+1)/hp->nt);
1257	–	if (j > 0) {
1258	–	d = dp[-1].r;
1259	–	if (dp[0].r > d) d = dp[0].r;
1260	–	mag1 += d * (nextsine - lastsine) *
1261	–	(b - bright(dp[-1].v));
1262	–	} else {
1263	–	d = dp[hp->np-1].r;
1264	–	if (dp[0].r > d) d = dp[0].r;
1265	–	mag1 += d * (nextsine - lastsine) *
1266	–	(b - bright(dp[hp->np-1].v));
1267	–	}
1268	–	dp += hp->np;
1269	–	lastsine = nextsine;
1270	–	}
1271	–	mag0 = 2.0PI / hp->np;
1272	–	phi = 2.0PI (double)j/hp->np;
1273	–	cosp = tcos(phi); sinp = tsin(phi);
1274	–	xd += mag0cosp - mag1sinp;
1275	–	yd += mag0sinp + mag1cosp;
1276	–	}
1277	–	for (i = 0; i < 3; i++)
1278	–	gv[i] = (xdhp->ux[i] + ydhp->uy[i])(hp->nthp->np)/PI;
1279	–	}
1280	–
1281	–
1282	–	void
1283	–	dirgradient( /* compute direction gradient */
1284	–	FVECT gv,
1285	–	AMBSAMP da, / assumes standard ordering */
1286	–	AMBHEMI *hp
1287	–	)
1288	–	{
1289	–	int i, j;
1290	–	double mag;
1291	–	double phi, xd, yd;
1292	–	AMBSAMP *dp;
1293	–
1294	–	xd = yd = 0.0;
1295	–	for (j = 0; j < hp->np; j++) {
1296	–	dp = da + j;
1297	–	mag = 0.0;
1298	–	for (i = 0; i < hp->nt; i++) {
1299	–	#ifdef DEBUG
1300	–	if (dp->t != i \|\| dp->p != j)
1301	–	error(CONSISTENCY,
1302	–	"division order in dirgradient");
1303	–	#endif
1304	–	/* tan(t) */
1305	–	mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
1306	–	dp += hp->np;
1307	–	}
1308	–	phi = 2.0PI (j+.5)/hp->np + PI/2.0;
1309	–	xd += mag * tcos(phi);
1310	–	yd += mag * tsin(phi);
1311	–	}
1312	–	for (i = 0; i < 3; i++)
1313	–	gv[i] = xdhp->ux[i] + ydhp->uy[i];
1314	–	}
1315	–
1316	–	#endif /* ! NEWAMB */

Diff Legend

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

Comparing ray/src/rt/ambcomp.c (file contents): Revision 2.55 by greg, Fri May 9 16:05:09 2014 UTC vs. Revision 2.90 by greg, Wed Nov 15 18:02:52 2023 UTC

Diff Legend

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.55 by greg, Fri May 9 16:05:09 2014 UTC vs.
Revision 2.90 by greg, Wed Nov 15 18:02:52 2023 UTC