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.52 by greg, Wed May 7 21:45:13 2014 UTC vs.
Revision 2.105 by greg, Thu Aug 21 20:38:41 2025 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	>	#ifndef MINSDIST
28	>	#define MINSDIST        0.2     /* def. min. spacing = 1/5th division */
29	>	#endif
30
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	–
31		typedef struct {
52	–	COLOR   v;              /* hemisphere sample value */
53	–	float   d;              /* reciprocal distance (1/rt) */
32		FVECT   p;              /* intersection point */
33	+	float   d;              /* reciprocal distance */
34	+	SCOLOR  v;              /* hemisphere sample value */
35		} AMBSAMP;              /* sample value */
36
37		typedef struct {
38		RAY     *rp;            /* originating ray sample */
59	–	FVECT   ux, uy;         /* tangent axis unit vectors */
39		int     ns;             /* number of samples per axis */
40	<	COLOR   acoef;          /* division contribution coefficient */
40	>	int     sampOK;         /* acquired full sample set? */
41	>	int     atyp;           /* RAMBIENT or TAMBIENT */
42	>	SCOLOR  acoef;          /* division contribution coefficient */
43	>	SCOLOR  acol;           /* accumulated color */
44	>	FVECT   onrm;           /* oriented unperturbed surface normal */
45	>	FVECT   ux, uy;         /* tangent axis unit vectors */
46		AMBSAMP sa[1];          /* sample array (extends struct) */
47		}  AMBHEMI;             /* ambient sample hemisphere */
48
49	<	#define ambndx(h,i,j)   ((i)*(h)->ns + (j))
50	<	#define ambsam(h,i,j)   (h)->sa[ambndx(h,i,j)]
49	>	#define AI(h,i,j)       ((i)*(h)->ns + (j))
50	>	#define ambsam(h,i,j)   (h)->sa[AI(h,i,j)]
51
52		typedef struct {
53		FVECT   r_i, r_i1, e_i, rcp, rI2_eJ2;
54		double  I1, I2;
71	–	int     valid;
55		} FFTRI;                /* vectors and coefficients for Hessian calculation */
56
57
58	<	/* Get index for adjacent vertex */
58	>	#define XLOTSIZ         512             /* size of used car lot */
59	>	#define CFIRST          0               /* first corner */
60	>	#define COTHER          (CFIRST+4)      /* non-corner sample */
61	>	#define CMAXTARGET      (int)(XLOTSIZ*MINSDIST/(1-MINSDIST))
62	>
63		static int
64	<	adjacent_verti(AMBHEMI *hp, int i, int j, int dbit)
64	>	psample_class(double ss[2])             /* classify patch sample */
65		{
66	<	int     i0 = i*hp->ns + j;
67	<
68	<	switch (dbit) {
69	<	case VDB_y:     return(i0 - hp->ns);
70	<	case VDB_x:     return(i0 - 1);
71	<	case VDB_Xy:    return(i0 - hp->ns + 1);
72	<	case VDB_xY:    return(i0 + hp->ns - 1);
73	<	case VDB_X:     return(i0 + 1);
74	<	case VDB_Y:     return(i0 + hp->ns);
75	<	/* the following should never occur */
89	<	case VDB_xy:    return(i0 - hp->ns - 1);
90	<	case VDB_XY:    return(i0 + hp->ns + 1);
66	>	if (ss[0] < MINSDIST) {
67	>	if (ss[1] < MINSDIST)
68	>	return(CFIRST);
69	>	if (ss[1] > 1.-MINSDIST)
70	>	return(CFIRST+2);
71	>	} else if (ss[0] > 1.-MINSDIST) {
72	>	if (ss[1] < MINSDIST)
73	>	return(CFIRST+1);
74	>	if (ss[1] > 1.-MINSDIST)
75	>	return(CFIRST+3);
76		}
77	<	return(-1);
77	>	return(COTHER);                 /* not in a corner */
78		}
79
80	<
81	<	/* Get vertex direction bit for the opposite edge to complete triangle */
97	<	static int
98	<	vdb_edge(int db1, int db2)
80	>	static void
81	>	trade_patchsamp(double ss[2])           /* trade in problem patch position */
82		{
83	<	switch (db1) {
84	<	case VDB_x:     return(db2==VDB_y ? VDB_Xy : VDB_Y);
85	<	case VDB_y:     return(db2==VDB_x ? VDB_xY : VDB_X);
86	<	case VDB_X:     return(db2==VDB_Xy ? VDB_y : VDB_xY);
87	<	case VDB_Y:     return(db2==VDB_xY ? VDB_x : VDB_Xy);
88	<	case VDB_xY:    return(db2==VDB_x ? VDB_y : VDB_X);
89	<	case VDB_Xy:    return(db2==VDB_y ? VDB_x : VDB_Y);
83	>	static float    tradelot[XLOTSIZ][2];
84	>	static int      gterm[COTHER+1];
85	>	double          repl[2];
86	>	int             sclass, rclass;
87	>	int             x;
88	>	re_initialize:                          /* initialize lot? */
89	>	while (gterm[COTHER] < XLOTSIZ) {
90	>	tradelot[gterm[COTHER]][0] = frandom();
91	>	tradelot[gterm[COTHER]][1] = frandom();
92	>	++gterm[COTHER];
93		}
94	<	error(INTERNAL, "forbidden diagonal in vdb_edge()");
95	<	return(-1);
94	>	/* get trade-in candidate... */
95	>	sclass = psample_class(ss);     /* submitted corner or not? */
96	>	switch (sclass) {
97	>	case COTHER:                    /* trade mid-edge with corner/any */
98	>	x = irandom( gterm[COTHER-1] > CMAXTARGET
99	>	? gterm[COTHER-1] : XLOTSIZ );
100	>	break;
101	>	case CFIRST:                    /* kick out of first corner */
102	>	x = gterm[CFIRST] + irandom(XLOTSIZ - gterm[CFIRST]);
103	>	break;
104	>	default:                        /* kick out of 2nd-4th corner */
105	>	x = irandom(XLOTSIZ - (gterm[sclass] - gterm[sclass-1]));
106	>	x += (x >= gterm[sclass-1])*(gterm[sclass] - gterm[sclass-1]);
107	>	break;
108	>	}
109	>	if (x >= XLOTSIZ) {             /* uh-oh... trapped in a corner! */
110	>	memset(gterm, 0, sizeof(gterm));
111	>	goto re_initialize;
112	>	}
113	>	repl[0] = tradelot[x][0];       /* save selected replacement (result) */
114	>	repl[1] = tradelot[x][1];
115	>	/* identify replacement class */
116	>	for (rclass = CFIRST; rclass < COTHER; rclass++)
117	>	if (x < gterm[rclass])
118	>	break;          /* repark to keep classes grouped */
119	>	while (rclass > sclass) {       /* replacement group after submitted? */
120	>	tradelot[x][0] = tradelot[gterm[rclass-1]][0];
121	>	tradelot[x][1] = tradelot[gterm[rclass-1]][1];
122	>	x = gterm[--rclass]++;
123	>	}
124	>	while (rclass < sclass) {       /* replacement group before submitted? */
125	>	tradelot[x][0] = tradelot[--gterm[rclass]][0];
126	>	tradelot[x][1] = tradelot[gterm[rclass]][1];
127	>	x = gterm[rclass++];
128	>	}
129	>	tradelot[x][0] = ss[0];         /* complete the trade-in */
130	>	tradelot[x][1] = ss[1];
131	>	ss[0] = repl[0];
132	>	ss[1] = repl[1];
133		}
134
135	+	#undef XLOTSIZ
136	+	#undef COTHER
137	+	#undef CFIRST
138
139	<	static AMBHEMI *
140	<	inithemi(                       /* initialize sampling hemisphere */
141	<	COLOR   ac,
142	<	RAY     *r,
143	<	double  wt
139	>
140	>	static int
141	>	ambcollision(                           /* proposed direction collides? */
142	>	AMBHEMI *hp,
143	>	int     i,
144	>	int     j,
145	>	RREAL   spt[2]
146		)
147		{
148	<	AMBHEMI *hp;
149	<	double  d;
150	<	int     n, i;
151	<	/* set number of divisions */
152	<	if (ambacc <= FTINY &&
153	<	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
154	<	wt = d;                 /* avoid ray termination */
155	<	n = sqrt(ambdiv * wt) + 0.5;
156	<	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
157	<	if (n < i)
158	<	n = i;
159	<	/* allocate sampling array */
160	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
161	<	if (hp == NULL)
162	<	return(NULL);
163	<	hp->rp = r;
164	<	hp->ns = n;
165	<	/* assign coefficient */
166	<	copycolor(hp->acoef, ac);
167	<	d = 1.0/(n*n);
168	<	scalecolor(hp->acoef, d);
169	<	/* make tangent plane axes */
170	<	hp->uy[0] = 0.5 - frandom();
171	<	hp->uy[1] = 0.5 - frandom();
172	<	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);
148	>	int     ii, jj;
149	>	/* check existing neighbors */
150	>	for (ii = i-1; ii <= i+1; ii++) {
151	>	if (ii < 0) continue;
152	>	if (ii >= hp->ns) break;
153	>	for (jj = j-1; jj <= j+1; jj++) {
154	>	AMBSAMP *ap;
155	>	FVECT   avec;
156	>	double  dx, dy;
157	>	if (jj < 0) continue;
158	>	if (jj >= hp->ns) break;
159	>	if ((ii==i) & (jj==j)) continue;
160	>	ap = &ambsam(hp,ii,jj);
161	>	if (ap->d <= .5/FHUGE)
162	>	continue;       /* no one home */
163	>	VSUB(avec, ap->p, hp->rp->rop);
164	>	normalize(avec);        /* use diskworld distance */
165	>	dx = DOT(avec, hp->ux) - spt[0];
166	>	dy = DOT(avec, hp->uy) - spt[1];
167	>	if ((dx*dx + dy*dy)*(hp->ns*hp->ns) <
168	>	PI*MINSDIST*MINSDIST)
169	>	return(1);      /* too close */
170	>	}
171	>	}
172	>	return(0);                      /* nothing to worry about */
173		}
174
175
159	–	/* Sample ambient division and apply weighting coefficient */
176		static int
177	<	getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
177	>	ambsample(                              /* initial ambient division sample */
178	>	AMBHEMI *hp,
179	>	int     i,
180	>	int     j,
181	>	int     n
182	>	)
183		{
184	+	int     trade_ok = (!n & (hp->ns >= 4))*21;
185	+	AMBSAMP *ap = &ambsam(hp,i,j);
186	+	RAY     ar;
187		int     hlist[3], ii;
188	<	double  spt[2], zd;
188	>	double  ss[2];
189	>	RREAL   spt[2];
190	>	double  zd;
191	>	/* generate hemispherical sample */
192		/* ambient coefficient for weight */
193		if (ambacc > FTINY)
194	<	setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL);
194	>	setscolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
195		else
196	<	copycolor(arp->rcoef, hp->acoef);
197	<	if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0)
196	>	copyscolor(ar.rcoef, hp->acoef);
197	>	if (rayorigin(&ar, hp->atyp, hp->rp, ar.rcoef) < 0)
198		return(0);
199		if (ambacc > FTINY) {
200	<	multcolor(arp->rcoef, hp->acoef);
201	<	scalecolor(arp->rcoef, 1./AVGREFL);
200	>	smultscolor(ar.rcoef, hp->acoef);
201	>	scalescolor(ar.rcoef, 1./AVGREFL);
202		}
203		hlist[0] = hp->rp->rno;
204	<	hlist[1] = j;
205	<	hlist[2] = i;
206	<	multisamp(spt, 2, urand(ilhash(hlist,3)+n));
207	<	if (!n) {                       /* avoid border samples for n==0 */
208	<	if ((spt[0] < 0.1) | (spt[0] >= 0.9))
209	<	spt[0] = 0.1 + 0.8*frandom();
210	<	if ((spt[1] < 0.1) | (spt[1] >= 0.9))
211	<	spt[1] = 0.1 + 0.8*frandom();
204	>	hlist[1] = AI(hp,i,j);
205	>	hlist[2] = samplendx;
206	>	multisamp(ss, 2, urand(ilhash(hlist,3)+n));
207	>	square2disk(spt, (j+ss[1])/hp->ns, (i+ss[0])/hp->ns);
208	>	/* avoid coincident samples? */
209	>	while (trade_ok-- && ambcollision(hp, i, j, spt)) {
210	>	if (trade_ok) {
211	>	trade_patchsamp(ss);
212	>	} else {                /* punting... */
213	>	ss[0] = MINSDIST + (1-2*MINSDIST)*frandom();
214	>	ss[1] = MINSDIST + (1-2*MINSDIST)*frandom();
215	>	}
216	>	square2disk(spt, (j+ss[1])/hp->ns, (i+ss[0])/hp->ns);
217		}
186	–	SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
218		zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
219		for (ii = 3; ii--; )
220	<	arp->rdir[ii] = spt[0]*hp->ux[ii] +
220	>	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
221		spt[1]*hp->uy[ii] +
222	<	zd*hp->rp->ron[ii];
223	<	checknorm(arp->rdir);
224	<	dimlist[ndims++] = ambndx(hp,i,j) + 90171;
225	<	rayvalue(arp);                  /* evaluate ray */
226	<	ndims--;                        /* apply coefficient */
227	<	multcolor(arp->rcol, arp->rcoef);
222	>	zd*hp->onrm[ii];
223	>	checknorm(ar.rdir);
224	>	dimlist[ndims_inc()] = AI(hp,i,j) + 90171;
225	>	rayvalue(&ar);                  /* evaluate ray */
226	>	dec_ndims();
227	>	zd = raydistance(&ar);
228	>	if (zd <= FTINY)
229	>	return(0);              /* should never happen */
230	>	smultscolor(ar.rcol, ar.rcoef); /* apply coefficient */
231	>	if (zd*ap->d < 1.0)             /* new/closer distance? */
232	>	ap->d = 1.0/zd;
233	>	if (!n) {                       /* record first vertex & value */
234	>	if (zd > 10.0*thescene.cusize + 1000.)
235	>	zd = 10.0*thescene.cusize + 1000.;
236	>	VSUM(ap->p, ar.rorg, ar.rdir, zd);
237	>	copyscolor(ap->v, ar.rcol);
238	>	} else {                        /* else update recorded value */
239	>	sopscolor(hp->acol, -=, ap->v);
240	>	zd = 1.0/(double)(n+1);
241	>	scalescolor(ar.rcol, zd);
242	>	zd *= (double)n;
243	>	scalescolor(ap->v, zd);
244	>	saddscolor(ap->v, ar.rcol);
245	>	}
246	>	saddscolor(hp->acol, ap->v);    /* add to our sum */
247		return(1);
248		}
249
250
251	<	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 */
251	>	/* Estimate variance based on ambient division differences */
252		static float *
253		getambdiffs(AMBHEMI *hp)
254		{
255	<	float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
255	>	const double    normf = 1./(pbright(hp->acoef) + FTINY);
256	>	float   *earr = (float *)calloc(2*hp->ns*hp->ns, sizeof(float));
257		float   *ep;
258		AMBSAMP *ap;
259	<	double  b, d2;
259	>	double  b, b1, d2;
260		int     i, j;
261
262		if (earr == NULL)               /* out of memory? */
263		return(NULL);
264	<	/* compute squared neighbor diffs */
265	<	for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
264	>	/* sum squared neighbor diffs */
265	>	ap = hp->sa;
266	>	ep = earr + hp->ns*hp->ns;      /* original estimates to scratch */
267	>	for (i = 0; i < hp->ns; i++)
268		for (j = 0; j < hp->ns; j++, ap++, ep++) {
269	<	b = bright(ap[0].v);
269	>	b = pbright(ap[0].v);
270		if (i) {                /* from above */
271	<	d2 = b - bright(ap[-hp->ns].v);
272	<	d2 *= d2;
271	>	b1 = pbright(ap[-hp->ns].v);
272	>	d2 = b - b1;
273	>	d2 *= d2*normf/(b + b1 + FTINY);
274		ep[0] += d2;
275		ep[-hp->ns] += d2;
276		}
277	<	if (j) {                /* from behind */
278	<	d2 = b - bright(ap[-1].v);
279	<	d2 *= d2;
280	<	ep[0] += d2;
281	<	ep[-1] += d2;
282	<	}
277	>	if (!j) continue;
278	>	/* from behind */
279	>	b1 = pbright(ap[-1].v);
280	>	d2 = b - b1;
281	>	d2 *= d2*normf/(b + b1 + FTINY);
282	>	ep[0] += d2;
283	>	ep[-1] += d2;
284	>	if (!i) continue;
285	>	/* diagonal */
286	>	b1 = pbright(ap[-hp->ns-1].v);
287	>	d2 = b - b1;
288	>	d2 *= d2*normf/(b + b1 + FTINY);
289	>	ep[0] += d2;
290	>	ep[-hp->ns-1] += d2;
291		}
292		/* correct for number of neighbors */
293	<	earr[0] *= 2.f;
294	<	earr[hp->ns-1] *= 2.f;
295	<	earr[(hp->ns-1)*hp->ns] *= 2.f;
296	<	earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f;
293	>	ep = earr + hp->ns*hp->ns;
294	>	ep[0] *= 6./3.;
295	>	ep[hp->ns-1] *= 6./3.;
296	>	ep[(hp->ns-1)*hp->ns] *= 6./3.;
297	>	ep[(hp->ns-1)*hp->ns + hp->ns-1] *= 6./3.;
298		for (i = 1; i < hp->ns-1; i++) {
299	<	earr[i*hp->ns] *= 4./3.;
300	<	earr[i*hp->ns + hp->ns-1] *= 4./3.;
299	>	ep[i*hp->ns] *= 6./5.;
300	>	ep[i*hp->ns + hp->ns-1] *= 6./5.;
301		}
302		for (j = 1; j < hp->ns-1; j++) {
303	<	earr[j] *= 4./3.;
304	<	earr[(hp->ns-1)*hp->ns + j] *= 4./3.;
303	>	ep[j] *= 6./5.;
304	>	ep[(hp->ns-1)*hp->ns + j] *= 6./5.;
305		}
306	+	/* blur final map to reduce bias */
307	+	for (i = 0; i < hp->ns-1; i++) {
308	+	float  *ep2;
309	+	ep = earr + i*hp->ns;
310	+	ep2 = ep + hp->ns*hp->ns;
311	+	for (j = 0; j < hp->ns-1; j++, ep++, ep2++) {
312	+	ep[0] += .5*ep2[0] + .125*(ep2[1] + ep2[hp->ns]);
313	+	ep[1] += .125*ep2[0];
314	+	ep[hp->ns] += .125*ep2[0];
315	+	}
316	+	}
317		return(earr);
318		}
319
320
321		/* Perform super-sampling on hemisphere (introduces bias) */
322		static void
323	<	ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
323	>	ambsupersamp(AMBHEMI *hp, int cnt)
324		{
325		float   *earr = getambdiffs(hp);
326	<	double  e2sum = 0.0;
276	<	AMBSAMP *ap;
277	<	RAY     ar;
278	<	double  asum[3];
326	>	double  e2rem = 0;
327		float   *ep;
328	<	int     i, j, n;
328	>	int     i, j, n, nss;
329
330		if (earr == NULL)               /* just skip calc. if no memory */
331		return;
332	<	/* add up estimated variances */
333	<	for (ep = earr + hp->ns*hp->ns; ep-- > earr; )
334	<	e2sum += *ep;
332	>	/* accumulate estimated variances */
333	>	for (ep = earr + hp->ns*hp->ns; ep > earr; )
334	>	e2rem += *--ep;
335		ep = earr;                      /* perform super-sampling */
336	<	for (ap = hp->sa, i = 0; i < hp->ns; i++)
337	<	for (j = 0; j < hp->ns; j++, ap++) {
338	<	int     nss = *ep/e2sum*cnt + frandom();
339	<	asum[0] = asum[1] = asum[2] = 0.0;
340	<	for (n = 1; n <= nss; n++) {
341	<	if (!getambsamp(&ar, hp, i, j, n)) {
342	<	nss = n-1;
343	<	break;
296	<	}
297	<	addcolor(asum, ar.rcol);
298	<	}
299	<	if (nss) {              /* update returned ambient value */
300	<	const double    ssf = 1./(nss + 1);
301	<	for (n = 3; n--; )
302	<	acol[n] += ssf*asum[n] +
303	<	(ssf - 1.)*colval(ap->v,n);
304	<	}
305	<	e2sum -= *ep++;         /* update remainders */
306	<	cnt -= nss;
336	>	for (i = 0; i < hp->ns; i++)
337	>	for (j = 0; j < hp->ns; j++) {
338	>	if (e2rem <= FTINY)
339	>	goto done;      /* nothing left to do */
340	>	nss = *ep/e2rem*cnt + frandom();
341	>	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
342	>	if (!--cnt) goto done;
343	>	e2rem -= *ep++;         /* update remainder */
344		}
345	+	done:
346		free(earr);
347		}
348
349
350	<	/* Compute vertex flags, indicating farthest in each direction */
351	<	static uby8 *
352	<	vertex_flags(AMBHEMI *hp)
350	>	static AMBHEMI *
351	>	samp_hemi(                              /* sample indirect hemisphere */
352	>	SCOLOR  rcol,
353	>	RAY     *r,
354	>	double  wt
355	>	)
356		{
357	<	uby8    *vflags = (uby8 *)calloc(hp->ns*hp->ns, sizeof(uby8));
358	<	uby8    *vf;
359	<	AMBSAMP *ap;
360	<	int     i, j;
357	>	int     backside = (wt < 0);
358	>	AMBHEMI *hp;
359	>	double  d;
360	>	int     n, i, j;
361	>	/* insignificance check */
362	>	d = sintens(rcol);
363	>	if (d <= FTINY)
364	>	return(NULL);
365	>	/* set number of divisions */
366	>	if (backside) wt = -wt;
367	>	if (ambacc <= FTINY &&
368	>	wt > (d *= 0.8*r->rweight/(ambdiv*minweight + 1e-20)))
369	>	wt = d;                 /* avoid ray termination */
370	>	n = sqrt(ambdiv * wt) + 0.5;
371	>	i = 1 + (MINADIV-1)*(ambacc > FTINY);
372	>	if (n < i)                      /* use minimum number of samples? */
373	>	n = i;
374	>	/* allocate sampling array */
375	>	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
376	>	if (hp == NULL)
377	>	error(SYSTEM, "out of memory in samp_hemi");
378
379	<	if (vflags == NULL)
380	<	error(SYSTEM, "out of memory in vertex_flags()");
381	<	vf = vflags;
382	<	ap = hp->sa;            /* compute farthest along first row */
383	<	for (j = 0; j < hp->ns-1; j++, vf++, ap++)
384	<	if (ap[0].d <= ap[1].d)
385	<	vf[0] |= 1<<VDB_X;
386	<	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;
379	>	if (backside) {
380	>	hp->atyp = TAMBIENT;
381	>	hp->onrm[0] = -r->ron[0];
382	>	hp->onrm[1] = -r->ron[1];
383	>	hp->onrm[2] = -r->ron[2];
384	>	} else {
385	>	hp->atyp = RAMBIENT;
386	>	VCOPY(hp->onrm, r->ron);
387		}
388	<	return(vflags);
388	>	hp->rp = r;
389	>	hp->ns = n;
390	>	scolorblack(hp->acol);
391	>	memset(hp->sa, 0, sizeof(AMBSAMP)*n*n);
392	>	hp->sampOK = 0;
393	>	/* assign coefficient */
394	>	copyscolor(hp->acoef, rcol);
395	>	d = 1.0/(n*n);
396	>	scalescolor(hp->acoef, d);
397	>	/* make tangent plane axes */
398	>	if (!getperpendicular(hp->ux, hp->onrm, 1))
399	>	error(CONSISTENCY, "bad ray direction in samp_hemi");
400	>	VCROSS(hp->uy, hp->onrm, hp->ux);
401	>	/* sample divisions */
402	>	for (i = hp->ns; i--; )
403	>	for (j = hp->ns; j--; )
404	>	hp->sampOK += ambsample(hp, i, j, 0);
405	>	copyscolor(rcol, hp->acol);
406	>	if (!hp->sampOK) {              /* utter failure? */
407	>	free(hp);
408	>	return(NULL);
409	>	}
410	>	if (hp->sampOK < hp->ns*hp->ns) {
411	>	hp->sampOK *= -1;       /* soft failure */
412	>	return(hp);
413	>	}
414	>	if (hp->sampOK <= MINADIV*MINADIV)
415	>	return(hp);             /* don't bother super-sampling */
416	>	n = ambssamp*wt + 0.5;
417	>	if (n >= 4*hp->ns) {            /* perform super-sampling? */
418	>	ambsupersamp(hp, n);
419	>	copyscolor(rcol, hp->acol);
420	>	}
421	>	return(hp);                     /* all is well */
422		}
423
424
425		/* Return brightness of farthest ambient sample */
426		static double
427	<	back_ambval(AMBHEMI *hp, int i, int j, int dbit1, int dbit2, const uby8 *vflags)
427	>	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
428		{
429	<	const int       v0 = ambndx(hp,i,j);
430	<	const int       tflags = (1<<dbit1 | 1<<dbit2);
431	<	int             v1, v2;
432	<
433	<	if ((vflags[v0] & tflags) == tflags)    /* is v0 the farthest? */
434	<	return(colval(hp->sa[v0].v,CIEY));
435	<	v1 = adjacent_verti(hp, i, j, dbit1);
436	<	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));
429	>	if (hp->sa[n1].d <= hp->sa[n2].d) {
430	>	if (hp->sa[n1].d <= hp->sa[n3].d)
431	>	return(hp->sa[n1].v[0]);
432	>	return(hp->sa[n3].v[0]);
433	>	}
434	>	if (hp->sa[n2].d <= hp->sa[n3].d)
435	>	return(hp->sa[n2].v[0]);
436	>	return(hp->sa[n3].v[0]);
437		}
438
439
440		/* Compute vectors and coefficients for Hessian/gradient calcs */
441		static void
442	<	comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags)
442	>	comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1)
443		{
444	<	const int       i0 = ambndx(hp,i,j);
445	<	double          rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
386	<	int             i1, ii;
444	>	double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
445	>	int     ii;
446
447	<	ftp->valid = 0;                 /* check if we can skip this edge */
448	<	ii = adjacent_trifl[dbit];
449	<	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);
447	>	VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop);
448	>	VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop);
449	>	VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p);
450		VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
451		rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
452		dot_e = DOT(ftp->e_i,ftp->e_i);
#	Line 410 | Line 460 | comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int
460		J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
461		for (ii = 3; ii--; )
462		ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
413	–	ftp->valid++;
463		}
464
465
#	Line 436 | Line 485 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
485		double  d1, d2, d3, d4;
486		double  I3, J3, K3;
487		int     i, j;
439	–
440	–	if (!ftp->valid) {              /* preemptive test */
441	–	memset(hess, 0, sizeof(FVECT)*3);
442	–	return;
443	–	}
488		/* compute intermediate coefficients */
489		d1 = 1.0/DOT(ftp->r_i,ftp->r_i);
490		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
#	Line 504 | Line 548 | comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
548		double  f1;
549		int     i;
550
507	–	if (!ftp->valid) {              /* preemptive test */
508	–	memset(grad, 0, sizeof(FVECT));
509	–	return;
510	–	}
551		f1 = 2.0*DOT(nrm, ftp->rcp);
552		VCROSS(ncp, nrm, ftp->e_i);
553		for (i = 3; i--; )
#	Line 537 | Line 577 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
577
578
579		/* Compute anisotropic radii and eigenvector directions */
580	<	static int
580	>	static void
581		eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
582		{
583		double  hess2[2][2];
#	Line 559 | Line 599 | eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
599		if (i == 1)                     /* double-root (circle) */
600		evalue[1] = evalue[0];
601		if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) |
602	<	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) )
603	<	error(INTERNAL, "bad eigenvalue calculation");
604	<
602	>	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) {
603	>	ra[0] = ra[1] = maxarad;
604	>	return;
605	>	}
606		if (evalue[0] > evalue[1]) {
607		ra[0] = sqrt(sqrt(4.0/evalue[0]));
608		ra[1] = sqrt(sqrt(4.0/evalue[1]));
#	Line 596 | Line 637 | ambHessian(                            /* anisotropic radii & pos. gradient */
637		static char     memerrmsg[] = "out of memory in ambHessian()";
638		FVECT           (*hessrow)[3] = NULL;
639		FVECT           *gradrow = NULL;
599	–	uby8            *vflags;
640		FVECT           hessian[3];
641		FVECT           gradient;
642		FFTRI           fftr;
#	Line 618 | Line 658 | ambHessian(                            /* anisotropic radii & pos. gradient */
658		error(SYSTEM, memerrmsg);
659		memset(gradient, 0, sizeof(gradient));
660		}
621	–	/* get vertex position flags */
622	–	vflags = vertex_flags(hp);
661		/* compute first row of edges */
662		for (j = 0; j < hp->ns-1; j++) {
663	<	comp_fftri(&fftr, hp, 0, j, VDB_X, vflags);
663	>	comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1));
664		if (hessrow != NULL)
665	<	comp_hessian(hessrow[j], &fftr, hp->rp->ron);
665	>	comp_hessian(hessrow[j], &fftr, hp->onrm);
666		if (gradrow != NULL)
667	<	comp_gradient(gradrow[j], &fftr, hp->rp->ron);
667	>	comp_gradient(gradrow[j], &fftr, hp->onrm);
668		}
669		/* sum each row of triangles */
670		for (i = 0; i < hp->ns-1; i++) {
671		FVECT       hesscol[3];     /* compute first vertical edge */
672		FVECT       gradcol;
673	<	comp_fftri(&fftr, hp, i, 0, VDB_Y, vflags);
673	>	comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0));
674		if (hessrow != NULL)
675	<	comp_hessian(hesscol, &fftr, hp->rp->ron);
675	>	comp_hessian(hesscol, &fftr, hp->onrm);
676		if (gradrow != NULL)
677	<	comp_gradient(gradcol, &fftr, hp->rp->ron);
677	>	comp_gradient(gradcol, &fftr, hp->onrm);
678		for (j = 0; j < hp->ns-1; j++) {
679		FVECT   hessdia[3];     /* compute triangle contributions */
680		FVECT   graddia;
681		double  backg;
682	<	backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags);
682	>	backg = back_ambval(hp, AI(hp,i,j),
683	>	AI(hp,i,j+1), AI(hp,i+1,j));
684		/* diagonal (inner) edge */
685	<	comp_fftri(&fftr, hp, i, j+1, VDB_xY, vflags);
685	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j));
686		if (hessrow != NULL) {
687	<	comp_hessian(hessdia, &fftr, hp->rp->ron);
687	>	comp_hessian(hessdia, &fftr, hp->onrm);
688		rev_hessian(hesscol);
689		add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
690		}
691		if (gradrow != NULL) {
692	<	comp_gradient(graddia, &fftr, hp->rp->ron);
692	>	comp_gradient(graddia, &fftr, hp->onrm);
693		rev_gradient(gradcol);
694		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
695		}
696		/* initialize edge in next row */
697	<	comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags);
697	>	comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j));
698		if (hessrow != NULL)
699	<	comp_hessian(hessrow[j], &fftr, hp->rp->ron);
699	>	comp_hessian(hessrow[j], &fftr, hp->onrm);
700		if (gradrow != NULL)
701	<	comp_gradient(gradrow[j], &fftr, hp->rp->ron);
701	>	comp_gradient(gradrow[j], &fftr, hp->onrm);
702		/* new column edge & paired triangle */
703	<	backg = back_ambval(hp, i+1, j+1, VDB_x, VDB_y, vflags);
704	<	comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags);
703	>	backg = back_ambval(hp, AI(hp,i+1,j+1),
704	>	AI(hp,i+1,j), AI(hp,i,j+1));
705	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1));
706		if (hessrow != NULL) {
707	<	comp_hessian(hesscol, &fftr, hp->rp->ron);
707	>	comp_hessian(hesscol, &fftr, hp->onrm);
708		rev_hessian(hessdia);
709		add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
710		if (i < hp->ns-2)
711		rev_hessian(hessrow[j]);
712		}
713		if (gradrow != NULL) {
714	<	comp_gradient(gradcol, &fftr, hp->rp->ron);
714	>	comp_gradient(gradcol, &fftr, hp->onrm);
715		rev_gradient(graddia);
716		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
717		if (i < hp->ns-2)
#	Line 682 | Line 722 | ambHessian(                            /* anisotropic radii & pos. gradient */
722		/* release row buffers */
723		if (hessrow != NULL) free(hessrow);
724		if (gradrow != NULL) free(gradrow);
685	–	free(vflags);
725
726		if (ra != NULL)                 /* extract eigenvectors & radii */
727		eigenvectors(uv, ra, hessian);
#	Line 708 | Line 747 | ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
747		/* use vector for azimuth + 90deg */
748		VSUB(vd, ap->p, hp->rp->rop);
749		/* brightness over cosine factor */
750	<	gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
750	>	gfact = ap->v[0] / DOT(hp->onrm, vd);
751		/* sine = proj_radius/vd_length */
752		dgsum[0] -= DOT(uv[1], vd) * gfact;
753		dgsum[1] += DOT(uv[0], vd) * gfact;
#	Line 723 | Line 762 | static uint32
762		ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
763		{
764		const double    max_d = 1.0/(minarad*ambacc + 0.001);
765	<	const double    ang_res = 0.5*PI/(hp->ns-1);
766	<	const double    ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY));
765	>	const double    ang_res = 0.5*PI/hp->ns;
766	>	const double    ang_step = ang_res/((int)(16/PI*ang_res) + 1.01);
767		double          avg_d = 0;
768		uint32          flgs = 0;
769	+	FVECT           vec;
770	+	double          u, v;
771	+	double          ang, a1;
772		int             i, j;
773		/* don't bother for a few samples */
774	<	if (hp->ns < 12)
774	>	if (hp->ns < 8)
775		return(0);
776		/* check distances overhead */
777		for (i = hp->ns*3/4; i-- > hp->ns>>2; )
#	Line 744 | Line 786 | ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, c
786		for (i = 0; i < hp->ns; i++)
787		for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) {
788		AMBSAMP *ap = &ambsam(hp,i,j);
747	–	FVECT   vec;
748	–	double  u, v;
749	–	double  ang, a1;
750	–	int     abp;
789		if ((ap->d <= FTINY) | (ap->d >= max_d))
790		continue;       /* too far or too near */
791		VSUB(vec, ap->p, hp->rp->rop);
792	<	u = DOT(vec, uv[0]) * ap->d;
793	<	v = DOT(vec, uv[1]) * ap->d;
794	<	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0)
792	>	u = DOT(vec, uv[0]);
793	>	v = DOT(vec, uv[1]);
794	>	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v)
795		continue;       /* occluder outside ellipse */
796		ang = atan2a(v, u);     /* else set direction flags */
797	<	for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step)
797	>	for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step)
798		flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0)));
799		}
800		return(flgs);
#	Line 765 | Line 803 | ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, c
803
804		int
805		doambient(                              /* compute ambient component */
806	<	COLOR   rcol,                   /* input/output color */
806	>	SCOLOR  rcol,                   /* input/output color */
807		RAY     *r,
808	<	double  wt,
808	>	double  wt,                     /* negative for back side */
809		FVECT   uv[2],                  /* returned (optional) */
810		float   ra[2],                  /* returned (optional) */
811		float   pg[2],                  /* returned (optional) */
#	Line 775 | Line 813 | doambient(                             /* compute ambient component */
813		uint32  *crlp                   /* returned (optional) */
814		)
815		{
816	<	AMBHEMI *hp = inithemi(rcol, r, wt);
779	<	int     cnt;
816	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
817		FVECT   my_uv[2];
818	<	double  d, K, acol[3];
818	>	double  d, K;
819		AMBSAMP *ap;
820	<	int     i, j;
821	<	/* check/initialize */
785	<	if (hp == NULL)
786	<	return(0);
820	>	int     i;
821	>	/* clear return values */
822		if (uv != NULL)
823		memset(uv, 0, sizeof(FVECT)*2);
824		if (ra != NULL)
#	Line 794 | Line 829 | doambient(                             /* compute ambient component */
829		dg[0] = dg[1] = 0.0;
830		if (crlp != NULL)
831		*crlp = 0;
832	<	/* sample the hemisphere */
833	<	acol[0] = acol[1] = acol[2] = 0.0;
834	<	cnt = 0;
835	<	for (i = hp->ns; i--; )
836	<	for (j = hp->ns; j--; )
837	<	if ((ap = ambsample(hp, i, j)) != NULL) {
838	<	addcolor(acol, ap->v);
804	<	++cnt;
805	<	}
806	<	if (!cnt) {
807	<	setcolor(rcol, 0.0, 0.0, 0.0);
808	<	free(hp);
809	<	return(0);              /* no valid samples */
832	>	if (hp == NULL)                 /* sampling falure? */
833	>	return(0);
834	>
835	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) ||
836	>	(hp->sampOK < 0) | (hp->ns < MINADIV)) {
837	>	free(hp);               /* Hessian not requested/possible */
838	>	return(-1);             /* value-only return value */
839		}
840	<	if (cnt < hp->ns*hp->ns) {      /* incomplete sampling? */
841	<	copycolor(rcol, acol);
813	<	free(hp);
814	<	return(-1);             /* return value w/o Hessian */
815	<	}
816	<	cnt = ambssamp*wt + 0.5;        /* perform super-sampling? */
817	<	if (cnt > 8)
818	<	ambsupersamp(acol, hp, cnt);
819	<	copycolor(rcol, acol);          /* final indirect irradiance/PI */
820	<	if ((ra == NULL) & (pg == NULL) & (dg == NULL)) {
821	<	free(hp);
822	<	return(-1);             /* no radius or gradient calc. */
823	<	}
824	<	if ((d = bright(acol)) > FTINY) {       /* normalize Y values */
825	<	d = 0.99*(hp->ns*hp->ns)/d;
840	>	if ((d = scolor_mean(rcol)) > FTINY) {
841	>	d = 0.99*(hp->ns*hp->ns)/d;     /* normalize avg. values */
842		K = 0.01;
843		} else {                        /* or fall back on geometric Hessian */
844		K = 1.0;
845		pg = NULL;
846		dg = NULL;
847	+	crlp = NULL;
848		}
849	<	ap = hp->sa;                    /* relative Y channel from here on... */
849	>	ap = hp->sa;                    /* single channel from here on... */
850		for (i = hp->ns*hp->ns; i--; ap++)
851	<	colval(ap->v,CIEY) = bright(ap->v)*d + K;
851	>	ap->v[0] = scolor_mean(ap->v)*d + K;
852
853		if (uv == NULL)                 /* make sure we have axis pointers */
854		uv = my_uv;
#	Line 855 | Line 872 | doambient(                             /* compute ambient component */
872		if (ra[1] < minarad)
873		ra[1] = minarad;
874		}
875	<	ra[0] *= d = 1.0/sqrt(sqrt(wt));
875	>	ra[0] *= d = 1.0/sqrt(fabs(wt));
876		if ((ra[1] *= d) > 2.0*ra[0])
877		ra[1] = 2.0*ra[0];
878		if (ra[1] > maxarad) {
#	Line 863 | Line 880 | doambient(                             /* compute ambient component */
880		if (ra[0] > maxarad)
881		ra[0] = maxarad;
882		}
883	<	if (crlp != NULL)       /* flag encroached directions */
883	>	/* flag encroached directions */
884	>	if (crlp != NULL)       /* XXX doesn't update with changes to ambacc */
885		*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc);
886		if (pg != NULL) {       /* cap gradient if necessary */
887		d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];
#	Line 877 | Line 895 | doambient(                             /* compute ambient component */
895		free(hp);                       /* clean up and return */
896		return(1);
897		}
880	–
881	–
882	–	#else /* ! NEWAMB */
883	–
884	–
885	–	void
886	–	inithemi(                       /* initialize sampling hemisphere */
887	–	AMBHEMI  *hp,
888	–	COLOR ac,
889	–	RAY  *r,
890	–	double  wt
891	–	)
892	–	{
893	–	double  d;
894	–	int  i;
895	–	/* set number of divisions */
896	–	if (ambacc <= FTINY &&
897	–	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
898	–	wt = d;                 /* avoid ray termination */
899	–	hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
900	–	i = ambacc > FTINY ? 3 : 1;     /* minimum number of samples */
901	–	if (hp->nt < i)
902	–	hp->nt = i;
903	–	hp->np = PI * hp->nt + 0.5;
904	–	/* set number of super-samples */
905	–	hp->ns = ambssamp * wt + 0.5;
906	–	/* assign coefficient */
907	–	copycolor(hp->acoef, ac);
908	–	d = 1.0/(hp->nt*hp->np);
909	–	scalecolor(hp->acoef, d);
910	–	/* make axes */
911	–	VCOPY(hp->uz, r->ron);
912	–	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
913	–	for (i = 0; i < 3; i++)
914	–	if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
915	–	break;
916	–	if (i >= 3)
917	–	error(CONSISTENCY, "bad ray direction in inithemi");
918	–	hp->uy[i] = 1.0;
919	–	fcross(hp->ux, hp->uy, hp->uz);
920	–	normalize(hp->ux);
921	–	fcross(hp->uy, hp->uz, hp->ux);
922	–	}
923	–
924	–
925	–	int
926	–	divsample(                              /* sample a division */
927	–	AMBSAMP  *dp,
928	–	AMBHEMI  *h,
929	–	RAY  *r
930	–	)
931	–	{
932	–	RAY  ar;
933	–	int  hlist[3];
934	–	double  spt[2];
935	–	double  xd, yd, zd;
936	–	double  b2;
937	–	double  phi;
938	–	int  i;
939	–	/* ambient coefficient for weight */
940	–	if (ambacc > FTINY)
941	–	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
942	–	else
943	–	copycolor(ar.rcoef, h->acoef);
944	–	if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
945	–	return(-1);
946	–	if (ambacc > FTINY) {
947	–	multcolor(ar.rcoef, h->acoef);
948	–	scalecolor(ar.rcoef, 1./AVGREFL);
949	–	}
950	–	hlist[0] = r->rno;
951	–	hlist[1] = dp->t;
952	–	hlist[2] = dp->p;
953	–	multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
954	–	zd = sqrt((dp->t + spt[0])/h->nt);
955	–	phi = 2.0*PI * (dp->p + spt[1])/h->np;
956	–	xd = tcos(phi) * zd;
957	–	yd = tsin(phi) * zd;
958	–	zd = sqrt(1.0 - zd*zd);
959	–	for (i = 0; i < 3; i++)
960	–	ar.rdir[i] =    xd*h->ux[i] +
961	–	yd*h->uy[i] +
962	–	zd*h->uz[i];
963	–	checknorm(ar.rdir);
964	–	dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
965	–	rayvalue(&ar);
966	–	ndims--;
967	–	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
968	–	addcolor(dp->v, ar.rcol);
969	–	/* use rt to improve gradient calc */
970	–	if (ar.rt > FTINY && ar.rt < FHUGE)
971	–	dp->r += 1.0/ar.rt;
972	–	/* (re)initialize error */
973	–	if (dp->n++) {
974	–	b2 = bright(dp->v)/dp->n - bright(ar.rcol);
975	–	b2 = b2*b2 + dp->k*((dp->n-1)*(dp->n-1));
976	–	dp->k = b2/(dp->n*dp->n);
977	–	} else
978	–	dp->k = 0.0;
979	–	return(0);
980	–	}
981	–
982	–
983	–	static int
984	–	ambcmp(                                 /* decreasing order */
985	–	const void *p1,
986	–	const void *p2
987	–	)
988	–	{
989	–	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
990	–	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
991	–
992	–	if (d1->k < d2->k)
993	–	return(1);
994	–	if (d1->k > d2->k)
995	–	return(-1);
996	–	return(0);
997	–	}
998	–
999	–
1000	–	static int
1001	–	ambnorm(                                /* standard order */
1002	–	const void *p1,
1003	–	const void *p2
1004	–	)
1005	–	{
1006	–	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
1007	–	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
1008	–	int     c;
1009	–
1010	–	if ( (c = d1->t - d2->t) )
1011	–	return(c);
1012	–	return(d1->p - d2->p);
1013	–	}
1014	–
1015	–
1016	–	double
1017	–	doambient(                              /* compute ambient component */
1018	–	COLOR  rcol,
1019	–	RAY  *r,
1020	–	double  wt,
1021	–	FVECT  pg,
1022	–	FVECT  dg
1023	–	)
1024	–	{
1025	–	double  b, d=0;
1026	–	AMBHEMI  hemi;
1027	–	AMBSAMP  *div;
1028	–	AMBSAMP  dnew;
1029	–	double  acol[3];
1030	–	AMBSAMP  *dp;
1031	–	double  arad;
1032	–	int  divcnt;
1033	–	int  i, j;
1034	–	/* initialize hemisphere */
1035	–	inithemi(&hemi, rcol, r, wt);
1036	–	divcnt = hemi.nt * hemi.np;
1037	–	/* initialize */
1038	–	if (pg != NULL)
1039	–	pg[0] = pg[1] = pg[2] = 0.0;
1040	–	if (dg != NULL)
1041	–	dg[0] = dg[1] = dg[2] = 0.0;
1042	–	setcolor(rcol, 0.0, 0.0, 0.0);
1043	–	if (divcnt == 0)
1044	–	return(0.0);
1045	–	/* allocate super-samples */
1046	–	if (hemi.ns > 0 || pg != NULL || dg != NULL) {
1047	–	div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP));
1048	–	if (div == NULL)
1049	–	error(SYSTEM, "out of memory in doambient");
1050	–	} else
1051	–	div = NULL;
1052	–	/* sample the divisions */
1053	–	arad = 0.0;
1054	–	acol[0] = acol[1] = acol[2] = 0.0;
1055	–	if ((dp = div) == NULL)
1056	–	dp = &dnew;
1057	–	divcnt = 0;
1058	–	for (i = 0; i < hemi.nt; i++)
1059	–	for (j = 0; j < hemi.np; j++) {
1060	–	dp->t = i; dp->p = j;
1061	–	setcolor(dp->v, 0.0, 0.0, 0.0);
1062	–	dp->r = 0.0;
1063	–	dp->n = 0;
1064	–	if (divsample(dp, &hemi, r) < 0) {
1065	–	if (div != NULL)
1066	–	dp++;
1067	–	continue;
1068	–	}
1069	–	arad += dp->r;
1070	–	divcnt++;
1071	–	if (div != NULL)
1072	–	dp++;
1073	–	else
1074	–	addcolor(acol, dp->v);
1075	–	}
1076	–	if (!divcnt) {
1077	–	if (div != NULL)
1078	–	free((void *)div);
1079	–	return(0.0);            /* no samples taken */
1080	–	}
1081	–	if (divcnt < hemi.nt*hemi.np) {
1082	–	pg = dg = NULL;         /* incomplete sampling */
1083	–	hemi.ns = 0;
1084	–	} else if (arad > FTINY && divcnt/arad < minarad) {
1085	–	hemi.ns = 0;            /* close enough */
1086	–	} else if (hemi.ns > 0) {       /* else perform super-sampling? */
1087	–	comperrs(div, &hemi);                   /* compute errors */
1088	–	qsort(div, divcnt, sizeof(AMBSAMP), ambcmp);    /* sort divs */
1089	–	/* super-sample */
1090	–	for (i = hemi.ns; i > 0; i--) {
1091	–	dnew = *div;
1092	–	if (divsample(&dnew, &hemi, r) < 0) {
1093	–	dp++;
1094	–	continue;
1095	–	}
1096	–	dp = div;               /* reinsert */
1097	–	j = divcnt < i ? divcnt : i;
1098	–	while (--j > 0 && dnew.k < dp[1].k) {
1099	–	*dp = *(dp+1);
1100	–	dp++;
1101	–	}
1102	–	*dp = dnew;
1103	–	}
1104	–	if (pg != NULL || dg != NULL)   /* restore order */
1105	–	qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
1106	–	}
1107	–	/* compute returned values */
1108	–	if (div != NULL) {
1109	–	arad = 0.0;             /* note: divcnt may be < nt*np */
1110	–	for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
1111	–	arad += dp->r;
1112	–	if (dp->n > 1) {
1113	–	b = 1.0/dp->n;
1114	–	scalecolor(dp->v, b);
1115	–	dp->r *= b;
1116	–	dp->n = 1;
1117	–	}
1118	–	addcolor(acol, dp->v);
1119	–	}
1120	–	b = bright(acol);
1121	–	if (b > FTINY) {
1122	–	b = 1.0/b;      /* compute & normalize gradient(s) */
1123	–	if (pg != NULL) {
1124	–	posgradient(pg, div, &hemi);
1125	–	for (i = 0; i < 3; i++)
1126	–	pg[i] *= b;
1127	–	}
1128	–	if (dg != NULL) {
1129	–	dirgradient(dg, div, &hemi);
1130	–	for (i = 0; i < 3; i++)
1131	–	dg[i] *= b;
1132	–	}
1133	–	}
1134	–	free((void *)div);
1135	–	}
1136	–	copycolor(rcol, acol);
1137	–	if (arad <= FTINY)
1138	–	arad = maxarad;
1139	–	else
1140	–	arad = (divcnt+hemi.ns)/arad;
1141	–	if (pg != NULL) {               /* reduce radius if gradient large */
1142	–	d = DOT(pg,pg);
1143	–	if (d*arad*arad > 1.0)
1144	–	arad = 1.0/sqrt(d);
1145	–	}
1146	–	if (arad < minarad) {
1147	–	arad = minarad;
1148	–	if (pg != NULL && d*arad*arad > 1.0) {  /* cap gradient */
1149	–	d = 1.0/arad/sqrt(d);
1150	–	for (i = 0; i < 3; i++)
1151	–	pg[i] *= d;
1152	–	}
1153	–	}
1154	–	if ((arad /= sqrt(wt)) > maxarad)
1155	–	arad = maxarad;
1156	–	return(arad);
1157	–	}
1158	–
1159	–
1160	–	void
1161	–	comperrs(                       /* compute initial error estimates */
1162	–	AMBSAMP  *da,   /* assumes standard ordering */
1163	–	AMBHEMI  *hp
1164	–	)
1165	–	{
1166	–	double  b, b2;
1167	–	int  i, j;
1168	–	AMBSAMP  *dp;
1169	–	/* sum differences from neighbors */
1170	–	dp = da;
1171	–	for (i = 0; i < hp->nt; i++)
1172	–	for (j = 0; j < hp->np; j++) {
1173	–	#ifdef  DEBUG
1174	–	if (dp->t != i || dp->p != j)
1175	–	error(CONSISTENCY,
1176	–	"division order in comperrs");
1177	–	#endif
1178	–	b = bright(dp[0].v);
1179	–	if (i > 0) {            /* from above */
1180	–	b2 = bright(dp[-hp->np].v) - b;
1181	–	b2 *= b2 * 0.25;
1182	–	dp[0].k += b2;
1183	–	dp[-hp->np].k += b2;
1184	–	}
1185	–	if (j > 0) {            /* from behind */
1186	–	b2 = bright(dp[-1].v) - b;
1187	–	b2 *= b2 * 0.25;
1188	–	dp[0].k += b2;
1189	–	dp[-1].k += b2;
1190	–	} else {                /* around */
1191	–	b2 = bright(dp[hp->np-1].v) - b;
1192	–	b2 *= b2 * 0.25;
1193	–	dp[0].k += b2;
1194	–	dp[hp->np-1].k += b2;
1195	–	}
1196	–	dp++;
1197	–	}
1198	–	/* divide by number of neighbors */
1199	–	dp = da;
1200	–	for (j = 0; j < hp->np; j++)            /* top row */
1201	–	(dp++)->k *= 1.0/3.0;
1202	–	if (hp->nt < 2)
1203	–	return;
1204	–	for (i = 1; i < hp->nt-1; i++)          /* central region */
1205	–	for (j = 0; j < hp->np; j++)
1206	–	(dp++)->k *= 0.25;
1207	–	for (j = 0; j < hp->np; j++)            /* bottom row */
1208	–	(dp++)->k *= 1.0/3.0;
1209	–	}
1210	–
1211	–
1212	–	void
1213	–	posgradient(                                    /* compute position gradient */
1214	–	FVECT  gv,
1215	–	AMBSAMP  *da,                   /* assumes standard ordering */
1216	–	AMBHEMI  *hp
1217	–	)
1218	–	{
1219	–	int  i, j;
1220	–	double  nextsine, lastsine, b, d;
1221	–	double  mag0, mag1;
1222	–	double  phi, cosp, sinp, xd, yd;
1223	–	AMBSAMP  *dp;
1224	–
1225	–	xd = yd = 0.0;
1226	–	for (j = 0; j < hp->np; j++) {
1227	–	dp = da + j;
1228	–	mag0 = mag1 = 0.0;
1229	–	lastsine = 0.0;
1230	–	for (i = 0; i < hp->nt; i++) {
1231	–	#ifdef  DEBUG
1232	–	if (dp->t != i || dp->p != j)
1233	–	error(CONSISTENCY,
1234	–	"division order in posgradient");
1235	–	#endif
1236	–	b = bright(dp->v);
1237	–	if (i > 0) {
1238	–	d = dp[-hp->np].r;
1239	–	if (dp[0].r > d) d = dp[0].r;
1240	–	/* sin(t)*cos(t)^2 */
1241	–	d *= lastsine * (1.0 - (double)i/hp->nt);
1242	–	mag0 += d*(b - bright(dp[-hp->np].v));
1243	–	}
1244	–	nextsine = sqrt((double)(i+1)/hp->nt);
1245	–	if (j > 0) {
1246	–	d = dp[-1].r;
1247	–	if (dp[0].r > d) d = dp[0].r;
1248	–	mag1 += d * (nextsine - lastsine) *
1249	–	(b - bright(dp[-1].v));
1250	–	} else {
1251	–	d = dp[hp->np-1].r;
1252	–	if (dp[0].r > d) d = dp[0].r;
1253	–	mag1 += d * (nextsine - lastsine) *
1254	–	(b - bright(dp[hp->np-1].v));
1255	–	}
1256	–	dp += hp->np;
1257	–	lastsine = nextsine;
1258	–	}
1259	–	mag0 *= 2.0*PI / hp->np;
1260	–	phi = 2.0*PI * (double)j/hp->np;
1261	–	cosp = tcos(phi); sinp = tsin(phi);
1262	–	xd += mag0*cosp - mag1*sinp;
1263	–	yd += mag0*sinp + mag1*cosp;
1264	–	}
1265	–	for (i = 0; i < 3; i++)
1266	–	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI;
1267	–	}
1268	–
1269	–
1270	–	void
1271	–	dirgradient(                                    /* compute direction gradient */
1272	–	FVECT  gv,
1273	–	AMBSAMP  *da,                   /* assumes standard ordering */
1274	–	AMBHEMI  *hp
1275	–	)
1276	–	{
1277	–	int  i, j;
1278	–	double  mag;
1279	–	double  phi, xd, yd;
1280	–	AMBSAMP  *dp;
1281	–
1282	–	xd = yd = 0.0;
1283	–	for (j = 0; j < hp->np; j++) {
1284	–	dp = da + j;
1285	–	mag = 0.0;
1286	–	for (i = 0; i < hp->nt; i++) {
1287	–	#ifdef  DEBUG
1288	–	if (dp->t != i || dp->p != j)
1289	–	error(CONSISTENCY,
1290	–	"division order in dirgradient");
1291	–	#endif
1292	–	/* tan(t) */
1293	–	mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
1294	–	dp += hp->np;
1295	–	}
1296	–	phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
1297	–	xd += mag * tcos(phi);
1298	–	yd += mag * tsin(phi);
1299	–	}
1300	–	for (i = 0; i < 3; i++)
1301	–	gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
1302	–	}
1303	–
1304	–	#endif  /* ! NEWAMB */

Diff Legend

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