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.38 by greg, Sat Apr 26 15:54:17 2014 UTC vs.
Revision 2.105 by greg, Thu Aug 21 20:38:41 2025 UTC

#	Line 8 | Line 8 | static const char      RCSid[] = "$Id$";
8		*      for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz
9		*      from ACM SIGGRAPH Asia 2012 conference proceedings.
10		*
11	+	*  Added book-keeping optimization to avoid calculations that would
12	+	*      cancel due to traversal both directions on edges that are adjacent
13	+	*      to same-valued triangles.  This cuts about half of Hessian math.
14	+	*
15		*  Declarations of external symbols in ambient.h
16		*/
17
#	Line 17 | Line 21 | static const char      RCSid[] = "$Id$";
21		#include  "ambient.h"
22		#include  "random.h"
23
24	<	#ifdef NEWAMB
24	>	#ifndef 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
31	<	extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
31	>	typedef struct {
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 */
26	–	FVECT   ux, uy;         /* tangent axis unit vectors */
39		int     ns;             /* number of samples per axis */
40	<	COLOR   acoef;          /* division contribution coefficient */
41	<	struct s_ambsamp {
42	<	COLOR   v;              /* hemisphere sample value */
43	<	FVECT   p;              /* intersection point */
44	<	} sa[1];                /* sample array (extends struct) */
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 ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (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;
#	Line 40 | Line 55 | typedef struct {
55		} FFTRI;                /* vectors and coefficients for Hessian calculation */
56
57
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	+	psample_class(double ss[2])             /* classify patch sample */
65	+	{
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(COTHER);                 /* not in a corner */
78	+	}
79	+
80	+	static void
81	+	trade_patchsamp(double ss[2])           /* trade in problem patch position */
82	+	{
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	+	/* 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	+
140	+	static int
141	+	ambcollision(                           /* proposed direction collides? */
142	+	AMBHEMI *hp,
143	+	int     i,
144	+	int     j,
145	+	RREAL   spt[2]
146	+	)
147	+	{
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	+
176	+	static int
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  ss[2];
189	+	RREAL   spt[2];
190	+	double  zd;
191	+	/* generate hemispherical sample */
192	+	/* ambient coefficient for weight */
193	+	if (ambacc > FTINY)
194	+	setscolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
195	+	else
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	+	smultscolor(ar.rcoef, hp->acoef);
201	+	scalescolor(ar.rcoef, 1./AVGREFL);
202	+	}
203	+	hlist[0] = hp->rp->rno;
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	+	}
218	+	zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
219	+	for (ii = 3; ii--; )
220	+	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
221	+	spt[1]*hp->uy[ii] +
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	+	/* Estimate variance based on ambient division differences */
252	+	static float *
253	+	getambdiffs(AMBHEMI *hp)
254	+	{
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, b1, d2;
260	+	int     i, j;
261	+
262	+	if (earr == NULL)               /* out of memory? */
263	+	return(NULL);
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 = pbright(ap[0].v);
270	+	if (i) {                /* from above */
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) 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	+	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	+	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	+	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(AMBHEMI *hp, int cnt)
324	+	{
325	+	float   *earr = getambdiffs(hp);
326	+	double  e2rem = 0;
327	+	float   *ep;
328	+	int     i, j, n, nss;
329	+
330	+	if (earr == NULL)               /* just skip calc. if no memory */
331	+	return;
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 (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		static AMBHEMI *
351	<	inithemi(                       /* initialize sampling hemisphere */
352	<	COLOR   ac,
351	>	samp_hemi(                              /* sample indirect hemisphere */
352	>	SCOLOR  rcol,
353		RAY     *r,
354		double  wt
355		)
356		{
357	+	int     backside = (wt < 0);
358		AMBHEMI *hp;
359		double  d;
360	<	int     n, i;
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*intens(ac)*r->rweight/(ambdiv*minweight)))
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 + 5*(ambacc > FTINY);     /* minimum number of samples */
372	<	if (n < i)
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) +
63	<	sizeof(struct s_ambsamp)*(n*n - 1));
375	>	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
376		if (hp == NULL)
377	<	return(NULL);
377	>	error(SYSTEM, "out of memory in samp_hemi");
378	>
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		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	<	copycolor(hp->acoef, ac);
394	>	copyscolor(hp->acoef, rcol);
395		d = 1.0/(n*n);
396	<	scalecolor(hp->acoef, d);
396	>	scalescolor(hp->acoef, d);
397		/* make tangent plane axes */
398	<	hp->uy[0] = 0.5 - frandom();
399	<	hp->uy[1] = 0.5 - frandom();
400	<	hp->uy[2] = 0.5 - frandom();
401	<	for (i = 3; i--; )
402	<	if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
403	<	break;
404	<	if (i < 0)
405	<	error(CONSISTENCY, "bad ray direction in inithemi");
406	<	hp->uy[i] = 1.0;
407	<	VCROSS(hp->ux, hp->uy, r->ron);
408	<	normalize(hp->ux);
409	<	VCROSS(hp->uy, r->ron, hp->ux);
410	<	/* we're ready to sample */
411	<	return(hp);
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	<	static struct s_ambsamp *
426	<	ambsample(                              /* sample an ambient direction */
427	<	AMBHEMI *hp,
93	<	int     i,
94	<	int     j
95	<	)
425	>	/* Return brightness of farthest ambient sample */
426	>	static double
427	>	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
428		{
429	<	struct s_ambsamp        *ap = &ambsamp(hp,i,j);
430	<	RAY                     ar;
431	<	double                  spt[2], zd;
432	<	int                     ii;
101	<	/* ambient coefficient for weight */
102	<	if (ambacc > FTINY)
103	<	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
104	<	else
105	<	copycolor(ar.rcoef, hp->acoef);
106	<	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
107	<	goto badsample;
108	<	if (ambacc > FTINY) {
109	<	multcolor(ar.rcoef, hp->acoef);
110	<	scalecolor(ar.rcoef, 1./AVGREFL);
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	<	/* generate hemispherical sample */
435	<	SDsquare2disk(spt,      (i+.1+.8*frandom())/hp->ns,
436	<	(j+.1+.8*frandom())/hp->ns );
115	<	zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
116	<	for (ii = 3; ii--; )
117	<	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
118	<	spt[1]*hp->uy[ii] +
119	<	zd*hp->rp->ron[ii];
120	<	checknorm(ar.rdir);
121	<	dimlist[ndims++] = i*hp->ns + j + 90171;
122	<	rayvalue(&ar);                  /* evaluate ray */
123	<	ndims--;
124	<	/* limit vertex distance */
125	<	if (ar.rt > 10.0*thescene.cusize)
126	<	ar.rt = 10.0*thescene.cusize;
127	<	else if (ar.rt <= FTINY)        /* should never happen! */
128	<	goto badsample;
129	<	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
130	<	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
131	<	copycolor(ap->v, ar.rcol);
132	<	return(ap);
133	<	badsample:
134	<	setcolor(ap->v, 0., 0., 0.);
135	<	VCOPY(ap->p, hp->rp->rop);
136	<	return(NULL);
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, FVECT ap0, FVECT ap1, FVECT rop)
442	>	comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1)
443		{
444		double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
445	<	int     i;
445	>	int     ii;
446
447	<	VSUB(ftp->r_i, ap0, rop);
448	<	VSUB(ftp->r_i1, ap1, rop);
449	<	VSUB(ftp->e_i, ap1, ap0);
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 158 | Line 458 | comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop
458		ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
459		dot_e*ftp->I1 )*0.5*rdot_cp;
460		J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
461	<	for (i = 3; i--; )
462	<	ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i];
461	>	for (ii = 3; ii--; )
462	>	ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
463		}
464
465
#	Line 208 | Line 508 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
508		hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] +
509		2.0*J3*m4[i][j] );
510		hess[i][j] += d2*(i==j);
511	<	hess[i][j] *= 1.0/PI;
511	>	hess[i][j] *= -1.0/PI;
512		}
513		}
514
#	Line 230 | Line 530 | rev_hessian(FVECT hess[3])
530		/* Add to radiometric Hessian from the given triangle */
531		static void
532		add2hessian(FVECT hess[3], FVECT ehess1[3],
533	<	FVECT ehess2[3], FVECT ehess3[3], COLORV v)
533	>	FVECT ehess2[3], FVECT ehess3[3], double v)
534		{
535		int     i, j;
536
#	Line 251 | Line 551 | comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
551		f1 = 2.0*DOT(nrm, ftp->rcp);
552		VCROSS(ncp, nrm, ftp->e_i);
553		for (i = 3; i--; )
554	<	grad[i] = (-0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
554	>	grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
555		}
556
557
#	Line 267 | Line 567 | rev_gradient(FVECT grad)
567
568		/* Add to displacement gradient from the given triangle */
569		static void
570	<	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
570	>	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
571		{
572		int     i;
573
#	Line 276 | Line 576 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
576		}
577
578
279	–	/* Return brightness of furthest ambient sample */
280	–	static COLORV
281	–	back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2,
282	–	struct s_ambsamp *ap3, FVECT orig)
283	–	{
284	–	COLORV  vback;
285	–	FVECT   vec;
286	–	double  d2, d2best;
287	–
288	–	VSUB(vec, ap1->p, orig);
289	–	d2best = DOT(vec,vec);
290	–	vback = colval(ap1->v,CIEY);
291	–	VSUB(vec, ap2->p, orig);
292	–	d2 = DOT(vec,vec);
293	–	if (d2 > d2best) {
294	–	d2best = d2;
295	–	vback = colval(ap2->v,CIEY);
296	–	}
297	–	VSUB(vec, ap3->p, orig);
298	–	d2 = DOT(vec,vec);
299	–	if (d2 > d2best)
300	–	return(colval(ap3->v,CIEY));
301	–	return(vback);
302	–	}
303	–
304	–
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 325 | 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 385 | Line 660 | ambHessian(                            /* anisotropic radii & pos. gradient */
660		}
661		/* compute first row of edges */
662		for (j = 0; j < hp->ns-1; j++) {
663	<	comp_fftri(&fftr, ambsamp(hp,0,j).p,
389	<	ambsamp(hp,0,j+1).p, hp->rp->rop);
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, ambsamp(hp,i,0).p,
400	<	ambsamp(hp,i+1,0).p, hp->rp->rop);
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	<	COLORV  backg;
682	<	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
683	<	&ambsamp(hp,i+1,j), hp->rp->rop);
681	>	double  backg;
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, ambsamp(hp,i,j+1).p,
413	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
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 (gradient != NULL) {
692	<	comp_gradient(graddia, &fftr, hp->rp->ron);
691	>	if (gradrow != NULL) {
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, ambsamp(hp,i+1,j+1).p,
426	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
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(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
704	<	&ambsamp(hp,i+1,j), hp->rp->rop);
705	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
435	<	hp->rp->rop);
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 466 | Line 736 | ambHessian(                            /* anisotropic radii & pos. gradient */
736		static void
737		ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
738		{
739	<	struct s_ambsamp        *ap;
740	<	double                  dgsum[2];
741	<	int                     n;
742	<	FVECT                   vd;
743	<	double                  gfact;
739	>	AMBSAMP *ap;
740	>	double  dgsum[2];
741	>	int     n;
742	>	FVECT   vd;
743	>	double  gfact;
744
745		dgsum[0] = dgsum[1] = 0.0;      /* sum values times -tan(theta) */
746		for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
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);
751	<	/* -sine = -proj_radius/vd_length */
752	<	dgsum[0] += DOT(uv[1], vd) * gfact;
753	<	dgsum[1] -= DOT(uv[0], vd) * gfact;
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;
754		}
755		dg[0] = dgsum[0] / (hp->ns*hp->ns);
756		dg[1] = dgsum[1] / (hp->ns*hp->ns);
757		}
758
759
760	+	/* Compute potential light leak direction flags for cache value */
761	+	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;
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 < 8)
775	+	return(0);
776	+	/* check distances overhead */
777	+	for (i = hp->ns*3/4; i-- > hp->ns>>2; )
778	+	for (j = hp->ns*3/4; j-- > hp->ns>>2; )
779	+	avg_d += ambsam(hp,i,j).d;
780	+	avg_d *= 4.0/(hp->ns*hp->ns);
781	+	if (avg_d*r0 >= 1.0)            /* ceiling too low for corral? */
782	+	return(0);
783	+	if (avg_d >= max_d)             /* insurance */
784	+	return(0);
785	+	/* else circle around perimeter */
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);
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]);
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-ang_res; a1 <= ang+ang_res; a1 += ang_step)
798	+	flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0)));
799	+	}
800	+	return(flgs);
801	+	}
802	+
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) */
812	<	float   dg[2]                   /* returned (optional) */
812	>	float   dg[2],                  /* returned (optional) */
813	>	uint32  *crlp                   /* returned (optional) */
814		)
815		{
816	<	AMBHEMI                 *hp = inithemi(rcol, r, wt);
817	<	int                     cnt = 0;
818	<	FVECT                   my_uv[2];
819	<	double                  d, K, acol[3];
820	<	struct s_ambsamp        *ap;
821	<	int                     i, j;
507	<	/* check/initialize */
508	<	if (hp == NULL)
509	<	return(0);
816	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
817	>	FVECT   my_uv[2];
818	>	double  d, K;
819	>	AMBSAMP *ap;
820	>	int     i;
821	>	/* clear return values */
822		if (uv != NULL)
823		memset(uv, 0, sizeof(FVECT)*2);
824		if (ra != NULL)
#	Line 515 | Line 827 | doambient(                             /* compute ambient component */
827		pg[0] = pg[1] = 0.0;
828		if (dg != NULL)
829		dg[0] = dg[1] = 0.0;
830	<	/* sample the hemisphere */
831	<	acol[0] = acol[1] = acol[2] = 0.0;
832	<	for (i = hp->ns; i--; )
833	<	for (j = hp->ns; j--; )
834	<	if ((ap = ambsample(hp, i, j)) != NULL) {
835	<	addcolor(acol, ap->v);
836	<	++cnt;
837	<	}
838	<	if (!cnt) {
527	<	setcolor(rcol, 0.0, 0.0, 0.0);
528	<	free(hp);
529	<	return(0);              /* no valid samples */
830	>	if (crlp != NULL)
831	>	*crlp = 0;
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	<	copycolor(rcol, acol);          /* final indirect irradiance/PI */
841	<	if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
533	<	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
534	<	free(hp);
535	<	return(-1);             /* no radius or gradient calc. */
536	<	}
537	<	if (bright(acol) > FTINY) {     /* normalize Y values */
538	<	d = 0.99*cnt/bright(acol);
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 {                        /* geometric Hessian fall-back */
541	<	d = 0.0;
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 569 | 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 577 | Line 880 | doambient(                             /* compute ambient component */
880		if (ra[0] > maxarad)
881		ra[0] = maxarad;
882		}
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];
888		if (d > 1.0) {
#	Line 589 | Line 895 | doambient(                             /* compute ambient component */
895		free(hp);                       /* clean up and return */
896		return(1);
897		}
592	–
593	–
594	–	#else /* ! NEWAMB */
595	–
596	–
597	–	void
598	–	inithemi(                       /* initialize sampling hemisphere */
599	–	AMBHEMI  *hp,
600	–	COLOR ac,
601	–	RAY  *r,
602	–	double  wt
603	–	)
604	–	{
605	–	double  d;
606	–	int  i;
607	–	/* set number of divisions */
608	–	if (ambacc <= FTINY &&
609	–	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
610	–	wt = d;                 /* avoid ray termination */
611	–	hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
612	–	i = ambacc > FTINY ? 3 : 1;     /* minimum number of samples */
613	–	if (hp->nt < i)
614	–	hp->nt = i;
615	–	hp->np = PI * hp->nt + 0.5;
616	–	/* set number of super-samples */
617	–	hp->ns = ambssamp * wt + 0.5;
618	–	/* assign coefficient */
619	–	copycolor(hp->acoef, ac);
620	–	d = 1.0/(hp->nt*hp->np);
621	–	scalecolor(hp->acoef, d);
622	–	/* make axes */
623	–	VCOPY(hp->uz, r->ron);
624	–	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
625	–	for (i = 0; i < 3; i++)
626	–	if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
627	–	break;
628	–	if (i >= 3)
629	–	error(CONSISTENCY, "bad ray direction in inithemi");
630	–	hp->uy[i] = 1.0;
631	–	fcross(hp->ux, hp->uy, hp->uz);
632	–	normalize(hp->ux);
633	–	fcross(hp->uy, hp->uz, hp->ux);
634	–	}
635	–
636	–
637	–	int
638	–	divsample(                              /* sample a division */
639	–	AMBSAMP  *dp,
640	–	AMBHEMI  *h,
641	–	RAY  *r
642	–	)
643	–	{
644	–	RAY  ar;
645	–	int  hlist[3];
646	–	double  spt[2];
647	–	double  xd, yd, zd;
648	–	double  b2;
649	–	double  phi;
650	–	int  i;
651	–	/* ambient coefficient for weight */
652	–	if (ambacc > FTINY)
653	–	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
654	–	else
655	–	copycolor(ar.rcoef, h->acoef);
656	–	if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
657	–	return(-1);
658	–	if (ambacc > FTINY) {
659	–	multcolor(ar.rcoef, h->acoef);
660	–	scalecolor(ar.rcoef, 1./AVGREFL);
661	–	}
662	–	hlist[0] = r->rno;
663	–	hlist[1] = dp->t;
664	–	hlist[2] = dp->p;
665	–	multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
666	–	zd = sqrt((dp->t + spt[0])/h->nt);
667	–	phi = 2.0*PI * (dp->p + spt[1])/h->np;
668	–	xd = tcos(phi) * zd;
669	–	yd = tsin(phi) * zd;
670	–	zd = sqrt(1.0 - zd*zd);
671	–	for (i = 0; i < 3; i++)
672	–	ar.rdir[i] =    xd*h->ux[i] +
673	–	yd*h->uy[i] +
674	–	zd*h->uz[i];
675	–	checknorm(ar.rdir);
676	–	dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
677	–	rayvalue(&ar);
678	–	ndims--;
679	–	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
680	–	addcolor(dp->v, ar.rcol);
681	–	/* use rt to improve gradient calc */
682	–	if (ar.rt > FTINY && ar.rt < FHUGE)
683	–	dp->r += 1.0/ar.rt;
684	–	/* (re)initialize error */
685	–	if (dp->n++) {
686	–	b2 = bright(dp->v)/dp->n - bright(ar.rcol);
687	–	b2 = b2*b2 + dp->k*((dp->n-1)*(dp->n-1));
688	–	dp->k = b2/(dp->n*dp->n);
689	–	} else
690	–	dp->k = 0.0;
691	–	return(0);
692	–	}
693	–
694	–
695	–	static int
696	–	ambcmp(                                 /* decreasing order */
697	–	const void *p1,
698	–	const void *p2
699	–	)
700	–	{
701	–	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
702	–	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
703	–
704	–	if (d1->k < d2->k)
705	–	return(1);
706	–	if (d1->k > d2->k)
707	–	return(-1);
708	–	return(0);
709	–	}
710	–
711	–
712	–	static int
713	–	ambnorm(                                /* standard order */
714	–	const void *p1,
715	–	const void *p2
716	–	)
717	–	{
718	–	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
719	–	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
720	–	int     c;
721	–
722	–	if ( (c = d1->t - d2->t) )
723	–	return(c);
724	–	return(d1->p - d2->p);
725	–	}
726	–
727	–
728	–	double
729	–	doambient(                              /* compute ambient component */
730	–	COLOR  rcol,
731	–	RAY  *r,
732	–	double  wt,
733	–	FVECT  pg,
734	–	FVECT  dg
735	–	)
736	–	{
737	–	double  b, d=0;
738	–	AMBHEMI  hemi;
739	–	AMBSAMP  *div;
740	–	AMBSAMP  dnew;
741	–	double  acol[3];
742	–	AMBSAMP  *dp;
743	–	double  arad;
744	–	int  divcnt;
745	–	int  i, j;
746	–	/* initialize hemisphere */
747	–	inithemi(&hemi, rcol, r, wt);
748	–	divcnt = hemi.nt * hemi.np;
749	–	/* initialize */
750	–	if (pg != NULL)
751	–	pg[0] = pg[1] = pg[2] = 0.0;
752	–	if (dg != NULL)
753	–	dg[0] = dg[1] = dg[2] = 0.0;
754	–	setcolor(rcol, 0.0, 0.0, 0.0);
755	–	if (divcnt == 0)
756	–	return(0.0);
757	–	/* allocate super-samples */
758	–	if (hemi.ns > 0 || pg != NULL || dg != NULL) {
759	–	div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP));
760	–	if (div == NULL)
761	–	error(SYSTEM, "out of memory in doambient");
762	–	} else
763	–	div = NULL;
764	–	/* sample the divisions */
765	–	arad = 0.0;
766	–	acol[0] = acol[1] = acol[2] = 0.0;
767	–	if ((dp = div) == NULL)
768	–	dp = &dnew;
769	–	divcnt = 0;
770	–	for (i = 0; i < hemi.nt; i++)
771	–	for (j = 0; j < hemi.np; j++) {
772	–	dp->t = i; dp->p = j;
773	–	setcolor(dp->v, 0.0, 0.0, 0.0);
774	–	dp->r = 0.0;
775	–	dp->n = 0;
776	–	if (divsample(dp, &hemi, r) < 0) {
777	–	if (div != NULL)
778	–	dp++;
779	–	continue;
780	–	}
781	–	arad += dp->r;
782	–	divcnt++;
783	–	if (div != NULL)
784	–	dp++;
785	–	else
786	–	addcolor(acol, dp->v);
787	–	}
788	–	if (!divcnt) {
789	–	if (div != NULL)
790	–	free((void *)div);
791	–	return(0.0);            /* no samples taken */
792	–	}
793	–	if (divcnt < hemi.nt*hemi.np) {
794	–	pg = dg = NULL;         /* incomplete sampling */
795	–	hemi.ns = 0;
796	–	} else if (arad > FTINY && divcnt/arad < minarad) {
797	–	hemi.ns = 0;            /* close enough */
798	–	} else if (hemi.ns > 0) {       /* else perform super-sampling? */
799	–	comperrs(div, &hemi);                   /* compute errors */
800	–	qsort(div, divcnt, sizeof(AMBSAMP), ambcmp);    /* sort divs */
801	–	/* super-sample */
802	–	for (i = hemi.ns; i > 0; i--) {
803	–	dnew = *div;
804	–	if (divsample(&dnew, &hemi, r) < 0) {
805	–	dp++;
806	–	continue;
807	–	}
808	–	dp = div;               /* reinsert */
809	–	j = divcnt < i ? divcnt : i;
810	–	while (--j > 0 && dnew.k < dp[1].k) {
811	–	*dp = *(dp+1);
812	–	dp++;
813	–	}
814	–	*dp = dnew;
815	–	}
816	–	if (pg != NULL || dg != NULL)   /* restore order */
817	–	qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
818	–	}
819	–	/* compute returned values */
820	–	if (div != NULL) {
821	–	arad = 0.0;             /* note: divcnt may be < nt*np */
822	–	for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
823	–	arad += dp->r;
824	–	if (dp->n > 1) {
825	–	b = 1.0/dp->n;
826	–	scalecolor(dp->v, b);
827	–	dp->r *= b;
828	–	dp->n = 1;
829	–	}
830	–	addcolor(acol, dp->v);
831	–	}
832	–	b = bright(acol);
833	–	if (b > FTINY) {
834	–	b = 1.0/b;      /* compute & normalize gradient(s) */
835	–	if (pg != NULL) {
836	–	posgradient(pg, div, &hemi);
837	–	for (i = 0; i < 3; i++)
838	–	pg[i] *= b;
839	–	}
840	–	if (dg != NULL) {
841	–	dirgradient(dg, div, &hemi);
842	–	for (i = 0; i < 3; i++)
843	–	dg[i] *= b;
844	–	}
845	–	}
846	–	free((void *)div);
847	–	}
848	–	copycolor(rcol, acol);
849	–	if (arad <= FTINY)
850	–	arad = maxarad;
851	–	else
852	–	arad = (divcnt+hemi.ns)/arad;
853	–	if (pg != NULL) {               /* reduce radius if gradient large */
854	–	d = DOT(pg,pg);
855	–	if (d*arad*arad > 1.0)
856	–	arad = 1.0/sqrt(d);
857	–	}
858	–	if (arad < minarad) {
859	–	arad = minarad;
860	–	if (pg != NULL && d*arad*arad > 1.0) {  /* cap gradient */
861	–	d = 1.0/arad/sqrt(d);
862	–	for (i = 0; i < 3; i++)
863	–	pg[i] *= d;
864	–	}
865	–	}
866	–	if ((arad /= sqrt(wt)) > maxarad)
867	–	arad = maxarad;
868	–	return(arad);
869	–	}
870	–
871	–
872	–	void
873	–	comperrs(                       /* compute initial error estimates */
874	–	AMBSAMP  *da,   /* assumes standard ordering */
875	–	AMBHEMI  *hp
876	–	)
877	–	{
878	–	double  b, b2;
879	–	int  i, j;
880	–	AMBSAMP  *dp;
881	–	/* sum differences from neighbors */
882	–	dp = da;
883	–	for (i = 0; i < hp->nt; i++)
884	–	for (j = 0; j < hp->np; j++) {
885	–	#ifdef  DEBUG
886	–	if (dp->t != i || dp->p != j)
887	–	error(CONSISTENCY,
888	–	"division order in comperrs");
889	–	#endif
890	–	b = bright(dp[0].v);
891	–	if (i > 0) {            /* from above */
892	–	b2 = bright(dp[-hp->np].v) - b;
893	–	b2 *= b2 * 0.25;
894	–	dp[0].k += b2;
895	–	dp[-hp->np].k += b2;
896	–	}
897	–	if (j > 0) {            /* from behind */
898	–	b2 = bright(dp[-1].v) - b;
899	–	b2 *= b2 * 0.25;
900	–	dp[0].k += b2;
901	–	dp[-1].k += b2;
902	–	} else {                /* around */
903	–	b2 = bright(dp[hp->np-1].v) - b;
904	–	b2 *= b2 * 0.25;
905	–	dp[0].k += b2;
906	–	dp[hp->np-1].k += b2;
907	–	}
908	–	dp++;
909	–	}
910	–	/* divide by number of neighbors */
911	–	dp = da;
912	–	for (j = 0; j < hp->np; j++)            /* top row */
913	–	(dp++)->k *= 1.0/3.0;
914	–	if (hp->nt < 2)
915	–	return;
916	–	for (i = 1; i < hp->nt-1; i++)          /* central region */
917	–	for (j = 0; j < hp->np; j++)
918	–	(dp++)->k *= 0.25;
919	–	for (j = 0; j < hp->np; j++)            /* bottom row */
920	–	(dp++)->k *= 1.0/3.0;
921	–	}
922	–
923	–
924	–	void
925	–	posgradient(                                    /* compute position gradient */
926	–	FVECT  gv,
927	–	AMBSAMP  *da,                   /* assumes standard ordering */
928	–	AMBHEMI  *hp
929	–	)
930	–	{
931	–	int  i, j;
932	–	double  nextsine, lastsine, b, d;
933	–	double  mag0, mag1;
934	–	double  phi, cosp, sinp, xd, yd;
935	–	AMBSAMP  *dp;
936	–
937	–	xd = yd = 0.0;
938	–	for (j = 0; j < hp->np; j++) {
939	–	dp = da + j;
940	–	mag0 = mag1 = 0.0;
941	–	lastsine = 0.0;
942	–	for (i = 0; i < hp->nt; i++) {
943	–	#ifdef  DEBUG
944	–	if (dp->t != i || dp->p != j)
945	–	error(CONSISTENCY,
946	–	"division order in posgradient");
947	–	#endif
948	–	b = bright(dp->v);
949	–	if (i > 0) {
950	–	d = dp[-hp->np].r;
951	–	if (dp[0].r > d) d = dp[0].r;
952	–	/* sin(t)*cos(t)^2 */
953	–	d *= lastsine * (1.0 - (double)i/hp->nt);
954	–	mag0 += d*(b - bright(dp[-hp->np].v));
955	–	}
956	–	nextsine = sqrt((double)(i+1)/hp->nt);
957	–	if (j > 0) {
958	–	d = dp[-1].r;
959	–	if (dp[0].r > d) d = dp[0].r;
960	–	mag1 += d * (nextsine - lastsine) *
961	–	(b - bright(dp[-1].v));
962	–	} else {
963	–	d = dp[hp->np-1].r;
964	–	if (dp[0].r > d) d = dp[0].r;
965	–	mag1 += d * (nextsine - lastsine) *
966	–	(b - bright(dp[hp->np-1].v));
967	–	}
968	–	dp += hp->np;
969	–	lastsine = nextsine;
970	–	}
971	–	mag0 *= 2.0*PI / hp->np;
972	–	phi = 2.0*PI * (double)j/hp->np;
973	–	cosp = tcos(phi); sinp = tsin(phi);
974	–	xd += mag0*cosp - mag1*sinp;
975	–	yd += mag0*sinp + mag1*cosp;
976	–	}
977	–	for (i = 0; i < 3; i++)
978	–	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI;
979	–	}
980	–
981	–
982	–	void
983	–	dirgradient(                                    /* compute direction gradient */
984	–	FVECT  gv,
985	–	AMBSAMP  *da,                   /* assumes standard ordering */
986	–	AMBHEMI  *hp
987	–	)
988	–	{
989	–	int  i, j;
990	–	double  mag;
991	–	double  phi, xd, yd;
992	–	AMBSAMP  *dp;
993	–
994	–	xd = yd = 0.0;
995	–	for (j = 0; j < hp->np; j++) {
996	–	dp = da + j;
997	–	mag = 0.0;
998	–	for (i = 0; i < hp->nt; i++) {
999	–	#ifdef  DEBUG
1000	–	if (dp->t != i || dp->p != j)
1001	–	error(CONSISTENCY,
1002	–	"division order in dirgradient");
1003	–	#endif
1004	–	/* tan(t) */
1005	–	mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
1006	–	dp += hp->np;
1007	–	}
1008	–	phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
1009	–	xd += mag * tcos(phi);
1010	–	yd += mag * tsin(phi);
1011	–	}
1012	–	for (i = 0; i < 3; i++)
1013	–	gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
1014	–	}
1015	–
1016	–	#endif  /* ! NEWAMB */

Diff Legend

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