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 1.1 by greg, Fri Jun 7 10:03:52 1991 UTC vs.
Revision 2.26 by greg, Wed Apr 16 20:32:00 2014 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1991 Regents of the University of California */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ LBL";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		* Routines to compute "ambient" values using Monte Carlo
6	+	*
7	+	*  Declarations of external symbols in ambient.h
8		*/
9
10	<	#include  "ray.h"
10	>	#include "copyright.h"
11
12	+	#include  "ray.h"
13		#include  "ambient.h"
14	–
14		#include  "random.h"
15
16	<	typedef struct {
18	<	short  t, p;            /* theta, phi indices */
19	<	COLOR  v;               /* value sum */
20	<	float  k;               /* error contribution for this division */
21	<	int  n;                 /* number of subsamples */
22	<	}  AMBSAMP;             /* ambient division sample */
16	>	#ifdef NEWAMB
17
18	+	extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
19	+
20		typedef struct {
21	<	FVECT  ux, uy, uz;      /* x, y and z axis directions */
22	<	short  nt, np;          /* number of theta and phi directions */
21	>	RAY     *rp;            /* originating ray sample */
22	>	FVECT   ux, uy;         /* tangent axis directions */
23	>	int     ns;             /* number of samples per axis */
24	>	COLOR   acoef;          /* division contribution coefficient */
25	>	struct s_ambsamp {
26	>	COLOR   v;              /* hemisphere sample value */
27	>	float   p[3];           /* intersection point */
28	>	} sa[1];                /* sample array (extends struct) */
29		}  AMBHEMI;             /* ambient sample hemisphere */
30
31	<	extern double  sin(), cos(), sqrt();
31	>	#define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
32
33
34	<	static int
35	<	ambcmp(d1, d2)                          /* decreasing order */
36	<	AMBSAMP  *d1, *d2;
34	>	static AMBHEMI *
35	>	inithemi(                       /* initialize sampling hemisphere */
36	>	COLOR   ac,
37	>	RAY     *r,
38	>	double  wt
39	>	)
40		{
41	<	if (d1->k < d2->k)
42	<	return(1);
43	<	if (d1->k > d2->k)
44	<	return(-1);
45	<	return(0);
41	>	AMBHEMI *hp;
42	>	double  d;
43	>	int     n, i;
44	>	/* set number of divisions */
45	>	if (ambacc <= FTINY &&
46	>	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
47	>	wt = d;                 /* avoid ray termination */
48	>	n = sqrt(ambdiv * wt) + 0.5;
49	>	i = 1 + 4*(ambacc > FTINY);     /* minimum number of samples */
50	>	if (n < i)
51	>	n = i;
52	>	/* allocate sampling array */
53	>	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
54	>	sizeof(struct s_ambsamp)*(n*n - 1));
55	>	if (hp == NULL)
56	>	return(NULL);
57	>	hp->rp = r;
58	>	hp->ns = n;
59	>	/* assign coefficient */
60	>	copycolor(hp->acoef, ac);
61	>	d = 1.0/(n*n);
62	>	scalecolor(hp->acoef, d);
63	>	/* make tangent axes */
64	>	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
65	>	for (i = 0; i < 3; i++)
66	>	if (r->rn[i] < 0.6 && r->rn[i] > -0.6)
67	>	break;
68	>	if (i >= 3)
69	>	error(CONSISTENCY, "bad ray direction in inithemi()");
70	>	hp->uy[i] = 1.0;
71	>	VCROSS(hp->ux, hp->uy, r->rn);
72	>	normalize(hp->ux);
73	>	VCROSS(hp->uy, r->rn, hp->ux);
74	>	/* we're ready to sample */
75	>	return(hp);
76		}
77
78
79		static int
80	<	ambnorm(d1, d2)                         /* standard order */
81	<	AMBSAMP  *d1, *d2;
80	>	ambsample(                              /* sample an ambient direction */
81	>	AMBHEMI *hp,
82	>	int     i,
83	>	int     j,
84	>	)
85		{
86	<	register int  c;
86	>	struct s_ambsamp        *ap = &ambsamp(hp,i,j);
87	>	RAY                     ar;
88	>	int                     hlist[3];
89	>	double                  spt[2], dz;
90	>	int                     ii;
91	>	/* ambient coefficient for weight */
92	>	if (ambacc > FTINY)
93	>	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
94	>	else
95	>	copycolor(ar.rcoef, hp->acoef);
96	>	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) {
97	>	setcolor(ap->v, 0., 0., 0.);
98	>	ap->r = 0.;
99	>	return(0);              /* no sample taken */
100	>	}
101	>	if (ambacc > FTINY) {
102	>	multcolor(ar.rcoef, hp->acoef);
103	>	scalecolor(ar.rcoef, 1./AVGREFL);
104	>	}
105	>	/* generate hemispherical sample */
106	>	SDsquare2disk(spt, (i+frandom())/hp->ns, (j+frandom())/hp->ns);
107	>	zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
108	>	for (ii = 3; ii--; )
109	>	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
110	>	spt[1]*hp->uy[ii] +
111	>	zd*hp->rp->ron[ii];
112	>	checknorm(ar.rdir);
113	>	dimlist[ndims++] = i*hp->ns + j + 90171;
114	>	rayvalue(&ar);                  /* evaluate ray */
115	>	ndims--;
116	>	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
117	>	copycolor(ap->v, ar.rcol);
118	>	if (ar.rt > 20.0*maxarad)       /* limit vertex distance */
119	>	ar.rt = 20.0*maxarad;
120	>	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
121	>	return(1);
122	>	}
123
124	<	if (c = d1->t - d2->t)
125	<	return(c);
126	<	return(d1->p - d2->p);
124	>
125	>	static void
126	>	ambHessian(                             /* anisotropic radii & pos. gradient */
127	>	AMBHEMI *hp,
128	>	FVECT   uv[2],                  /* returned */
129	>	float   ra[2],                  /* returned */
130	>	float   pg[2]                   /* returned */
131	>	)
132	>	{
133	>	if (ra != NULL) {               /* compute Hessian-derived radii */
134	>	} else {                        /* else copy original tangent axes */
135	>	VCOPY(uv[0], hp->ux);
136	>	VCOPY(uv[1], hp->uy);
137	>	}
138	>	if (pg == NULL)                 /* no position gradient requested? */
139	>	return;
140		}
141
142	+	int
143	+	doambient(                              /* compute ambient component */
144	+	COLOR   rcol,                   /* input/output color */
145	+	RAY     *r,
146	+	double  wt,
147	+	FVECT   uv[2],                  /* returned */
148	+	float   ra[2],                  /* returned */
149	+	float   pg[2],                  /* returned */
150	+	float   dg[2]                   /* returned */
151	+	)
152	+	{
153	+	int                     cnt = 0;
154	+	FVECT                   my_uv[2];
155	+	AMBHEMI                 *hp;
156	+	double                  d, acol[3];
157	+	struct s_ambsamp        *ap;
158	+	int                     i, j;
159	+	/* initialize */
160	+	if ((hp = inithemi(rcol, r, wt)) == NULL)
161	+	return(0);
162	+	if (uv != NULL)
163	+	memset(uv, 0, sizeof(FVECT)*2);
164	+	if (ra != NULL)
165	+	ra[0] = ra[1] = 0.0;
166	+	if (pg != NULL)
167	+	pg[0] = pg[1] = 0.0;
168	+	if (dg != NULL)
169	+	dg[0] = dg[1] = 0.0;
170	+	/* sample the hemisphere */
171	+	acol[0] = acol[1] = acol[2] = 0.0;
172	+	for (i = hemi.ns; i--; )
173	+	for (j = hemi.ns; j--; )
174	+	if (ambsample(hp, i, j)) {
175	+	ap = &ambsamp(hp,i,j);
176	+	addcolor(acol, ap->v);
177	+	++cnt;
178	+	}
179	+	if (!cnt) {
180	+	setcolor(rcol, 0.0, 0.0, 0.0);
181	+	free(hp);
182	+	return(0);              /* no valid samples */
183	+	}
184	+	d = 1.0 / cnt;                  /* final indirect irradiance/PI */
185	+	acol[0] *= d; acol[1] *= d; acol[2] *= d;
186	+	copycolor(rcol, acol);
187	+	if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
188	+	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
189	+	free(hp);
190	+	return(-1);             /* no radius or gradient calc. */
191	+	}
192	+	d = 0.01 * bright(rcol);        /* add in 1% before Hessian comp. */
193	+	if (d < FTINY) d = FTINY;
194	+	ap = hp->sa;                    /* using Y channel from here on... */
195	+	for (i = hp->ns*hp->ns; i--; ap++)
196	+	colval(ap->v,CIEY) = bright(ap->v) + d;
197
198	<	static double
199	<	divsample(dp, h, r)                     /* sample a division */
200	<	register AMBSAMP  *dp;
201	<	AMBHEMI  *h;
202	<	RAY  *r;
198	>	if (uv == NULL)                 /* make sure we have axis pointers */
199	>	uv = my_uv;
200	>	/* compute radii & pos. gradient */
201	>	ambHessian(hp, uv, ra, pg);
202	>	if (dg != NULL)                 /* compute direction gradient */
203	>	ambdirgrad(hp, uv, dg);
204	>	if (ra != NULL) {               /* adjust/clamp radii */
205	>	d = pow(wt, -0.25);
206	>	if ((ra[0] *= d) > maxarad)
207	>	ra[0] = maxarad;
208	>	if ((ra[1] *= d) > 2.0*ra[0])
209	>	ra[1] = 2.0*ra[0];
210	>	}
211	>	free(hp);                       /* clean up and return */
212	>	return(1);
213	>	}
214	>
215	>
216	>	#else /* ! NEWAMB */
217	>
218	>
219	>	void
220	>	inithemi(                       /* initialize sampling hemisphere */
221	>	AMBHEMI  *hp,
222	>	COLOR ac,
223	>	RAY  *r,
224	>	double  wt
225	>	)
226		{
227	+	double  d;
228	+	int  i;
229	+	/* set number of divisions */
230	+	if (ambacc <= FTINY &&
231	+	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
232	+	wt = d;                 /* avoid ray termination */
233	+	hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
234	+	i = ambacc > FTINY ? 3 : 1;     /* minimum number of samples */
235	+	if (hp->nt < i)
236	+	hp->nt = i;
237	+	hp->np = PI * hp->nt + 0.5;
238	+	/* set number of super-samples */
239	+	hp->ns = ambssamp * wt + 0.5;
240	+	/* assign coefficient */
241	+	copycolor(hp->acoef, ac);
242	+	d = 1.0/(hp->nt*hp->np);
243	+	scalecolor(hp->acoef, d);
244	+	/* make axes */
245	+	VCOPY(hp->uz, r->ron);
246	+	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
247	+	for (i = 0; i < 3; i++)
248	+	if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
249	+	break;
250	+	if (i >= 3)
251	+	error(CONSISTENCY, "bad ray direction in inithemi");
252	+	hp->uy[i] = 1.0;
253	+	fcross(hp->ux, hp->uy, hp->uz);
254	+	normalize(hp->ux);
255	+	fcross(hp->uy, hp->uz, hp->ux);
256	+	}
257	+
258	+
259	+	int
260	+	divsample(                              /* sample a division */
261	+	AMBSAMP  *dp,
262	+	AMBHEMI  *h,
263	+	RAY  *r
264	+	)
265	+	{
266		RAY  ar;
267	<	int  hlist[4];
267	>	int  hlist[3];
268	>	double  spt[2];
269		double  xd, yd, zd;
270		double  b2;
271		double  phi;
272	<	register int  k;
273	<
274	<	if (rayorigin(&ar, r, AMBIENT, 0.5) < 0)
275	<	return(0.0);
272	>	int  i;
273	>	/* ambient coefficient for weight */
274	>	if (ambacc > FTINY)
275	>	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
276	>	else
277	>	copycolor(ar.rcoef, h->acoef);
278	>	if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
279	>	return(-1);
280	>	if (ambacc > FTINY) {
281	>	multcolor(ar.rcoef, h->acoef);
282	>	scalecolor(ar.rcoef, 1./AVGREFL);
283	>	}
284		hlist[0] = r->rno;
285		hlist[1] = dp->t;
286		hlist[2] = dp->p;
287	<	hlist[3] = 0;
288	<	zd = sqrt((dp->t+urand(ilhash(hlist,4)+dp->n))/h->nt);
289	<	hlist[3] = 1;
290	<	phi = 2.0*PI * (dp->p+urand(ilhash(hlist,4)+dp->n))/h->np;
291	<	xd = cos(phi) * zd;
79	<	yd = sin(phi) * zd;
287	>	multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
288	>	zd = sqrt((dp->t + spt[0])/h->nt);
289	>	phi = 2.0*PI * (dp->p + spt[1])/h->np;
290	>	xd = tcos(phi) * zd;
291	>	yd = tsin(phi) * zd;
292		zd = sqrt(1.0 - zd*zd);
293	<	for (k = 0; k < 3; k++)
294	<	ar.rdir[k] =    xd*h->ux[k] +
295	<	yd*h->uy[k] +
296	<	zd*h->uz[k];
297	<	dimlist[ndims++] = dp->t*h->np + dp->p + 38813;
293	>	for (i = 0; i < 3; i++)
294	>	ar.rdir[i] =    xd*h->ux[i] +
295	>	yd*h->uy[i] +
296	>	zd*h->uz[i];
297	>	checknorm(ar.rdir);
298	>	dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
299		rayvalue(&ar);
300		ndims--;
301	+	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
302		addcolor(dp->v, ar.rcol);
303	+	/* use rt to improve gradient calc */
304	+	if (ar.rt > FTINY && ar.rt < FHUGE)
305	+	dp->r += 1.0/ar.rt;
306		/* (re)initialize error */
307		if (dp->n++) {
308		b2 = bright(dp->v)/dp->n - bright(ar.rcol);
#	Line 93 | Line 310 | RAY  *r;
310		dp->k = b2/(dp->n*dp->n);
311		} else
312		dp->k = 0.0;
313	<	return(ar.rot);
313	>	return(0);
314		}
315
316
317	+	static int
318	+	ambcmp(                                 /* decreasing order */
319	+	const void *p1,
320	+	const void *p2
321	+	)
322	+	{
323	+	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
324	+	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
325	+
326	+	if (d1->k < d2->k)
327	+	return(1);
328	+	if (d1->k > d2->k)
329	+	return(-1);
330	+	return(0);
331	+	}
332	+
333	+
334	+	static int
335	+	ambnorm(                                /* standard order */
336	+	const void *p1,
337	+	const void *p2
338	+	)
339	+	{
340	+	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
341	+	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
342	+	int     c;
343	+
344	+	if ( (c = d1->t - d2->t) )
345	+	return(c);
346	+	return(d1->p - d2->p);
347	+	}
348	+
349	+
350		double
351	<	doambient(acol, r, pg, dg)              /* compute ambient component */
352	<	COLOR  acol;
353	<	RAY  *r;
354	<	FVECT  pg, dg;
351	>	doambient(                              /* compute ambient component */
352	>	COLOR  rcol,
353	>	RAY  *r,
354	>	double  wt,
355	>	FVECT  pg,
356	>	FVECT  dg
357	>	)
358		{
359	<	double  b, d;
359	>	double  b, d=0;
360		AMBHEMI  hemi;
361		AMBSAMP  *div;
362		AMBSAMP  dnew;
363	<	register AMBSAMP  *dp;
363	>	double  acol[3];
364	>	AMBSAMP  *dp;
365		double  arad;
366	<	int  ndivs, ns;
367	<	register int  i, j;
114	<	/* initialize color */
115	<	setcolor(acol, 0.0, 0.0, 0.0);
366	>	int  divcnt;
367	>	int  i, j;
368		/* initialize hemisphere */
369	<	inithemi(&hemi, r);
370	<	ndivs = hemi.nt * hemi.np;
371	<	if (ndivs == 0)
369	>	inithemi(&hemi, rcol, r, wt);
370	>	divcnt = hemi.nt * hemi.np;
371	>	/* initialize */
372	>	if (pg != NULL)
373	>	pg[0] = pg[1] = pg[2] = 0.0;
374	>	if (dg != NULL)
375	>	dg[0] = dg[1] = dg[2] = 0.0;
376	>	setcolor(rcol, 0.0, 0.0, 0.0);
377	>	if (divcnt == 0)
378		return(0.0);
379	<	/* set number of super-samples */
380	<	ns = ambssamp * r->rweight + 0.5;
381	<	if (ns > 0 || pg != NULL || dg != NULL) {
124	<	div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
379	>	/* allocate super-samples */
380	>	if (hemi.ns > 0 || pg != NULL || dg != NULL) {
381	>	div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP));
382		if (div == NULL)
383		error(SYSTEM, "out of memory in doambient");
384		} else
385		div = NULL;
386		/* sample the divisions */
387		arad = 0.0;
388	+	acol[0] = acol[1] = acol[2] = 0.0;
389		if ((dp = div) == NULL)
390		dp = &dnew;
391	+	divcnt = 0;
392		for (i = 0; i < hemi.nt; i++)
393		for (j = 0; j < hemi.np; j++) {
394		dp->t = i; dp->p = j;
395		setcolor(dp->v, 0.0, 0.0, 0.0);
396	+	dp->r = 0.0;
397		dp->n = 0;
398	<	if ((d = divsample(dp, &hemi, r)) == 0.0)
399	<	goto oopsy;
400	<	if (d < FHUGE)
401	<	arad += 1.0 / d;
398	>	if (divsample(dp, &hemi, r) < 0) {
399	>	if (div != NULL)
400	>	dp++;
401	>	continue;
402	>	}
403	>	arad += dp->r;
404	>	divcnt++;
405		if (div != NULL)
406		dp++;
407		else
408		addcolor(acol, dp->v);
409		}
410	<	if (ns > 0) {                   /* perform super-sampling */
411	<	comperrs(div, hemi);                    /* compute errors */
412	<	qsort(div, ndivs, sizeof(AMBSAMP), ambcmp);     /* sort divs */
413	<	dp = div + ndivs;                       /* skim excess */
414	<	for (i = ndivs; i > ns; i--) {
415	<	dp--;
416	<	addcolor(acol, dp->v);
417	<	}
410	>	if (!divcnt) {
411	>	if (div != NULL)
412	>	free((void *)div);
413	>	return(0.0);            /* no samples taken */
414	>	}
415	>	if (divcnt < hemi.nt*hemi.np) {
416	>	pg = dg = NULL;         /* incomplete sampling */
417	>	hemi.ns = 0;
418	>	} else if (arad > FTINY && divcnt/arad < minarad) {
419	>	hemi.ns = 0;            /* close enough */
420	>	} else if (hemi.ns > 0) {       /* else perform super-sampling? */
421	>	comperrs(div, &hemi);                   /* compute errors */
422	>	qsort(div, divcnt, sizeof(AMBSAMP), ambcmp);    /* sort divs */
423		/* super-sample */
424	<	for (i = ns; i > 0; i--) {
425	<	copystruct(&dnew, div);
426	<	if ((d = divsample(&dnew, &hemi)) == 0.0)
427	<	goto oopsy;
428	<	if (d < FHUGE)
429	<	arad += 1.0 / d;
430	<	/* reinsert */
431	<	dp = div;
164	<	j = ndivs < i ? ndivs : i;
424	>	for (i = hemi.ns; i > 0; i--) {
425	>	dnew = *div;
426	>	if (divsample(&dnew, &hemi, r) < 0) {
427	>	dp++;
428	>	continue;
429	>	}
430	>	dp = div;               /* reinsert */
431	>	j = divcnt < i ? divcnt : i;
432		while (--j > 0 && dnew.k < dp[1].k) {
433	<	copystruct(dp, dp+1);
433	>	*dp = *(dp+1);
434		dp++;
435		}
436	<	copystruct(dp, &dnew);
170	<	/* extract darkest */
171	<	if (i <= ndivs) {
172	<	dp = div + i-1;
173	<	if (dp->n > 1) {
174	<	b = 1.0/dp->n;
175	<	scalecolor(dp->v, b);
176	<	dp->n = 1;
177	<	}
178	<	addcolor(acol, dp->v);
179	<	}
436	>	*dp = dnew;
437		}
438	<	if (pg != NULL || dg != NULL)   /* reorder */
439	<	qsort(div, ndivs, sizeof(AMBSAMP), ambnorm);
438	>	if (pg != NULL || dg != NULL)   /* restore order */
439	>	qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
440		}
441		/* compute returned values */
442	<	if (pg != NULL)
443	<	posgradient(pg, div, &hemi);
444	<	if (dg != NULL)
445	<	dirgradient(dg, div, &hemi);
446	<	if (div != NULL)
447	<	free((char *)div);
448	<	b = 1.0/ndivs;
449	<	scalecolor(acol, b);
442	>	if (div != NULL) {
443	>	arad = 0.0;             /* note: divcnt may be < nt*np */
444	>	for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
445	>	arad += dp->r;
446	>	if (dp->n > 1) {
447	>	b = 1.0/dp->n;
448	>	scalecolor(dp->v, b);
449	>	dp->r *= b;
450	>	dp->n = 1;
451	>	}
452	>	addcolor(acol, dp->v);
453	>	}
454	>	b = bright(acol);
455	>	if (b > FTINY) {
456	>	b = 1.0/b;      /* compute & normalize gradient(s) */
457	>	if (pg != NULL) {
458	>	posgradient(pg, div, &hemi);
459	>	for (i = 0; i < 3; i++)
460	>	pg[i] *= b;
461	>	}
462	>	if (dg != NULL) {
463	>	dirgradient(dg, div, &hemi);
464	>	for (i = 0; i < 3; i++)
465	>	dg[i] *= b;
466	>	}
467	>	}
468	>	free((void *)div);
469	>	}
470	>	copycolor(rcol, acol);
471		if (arad <= FTINY)
194	–	arad = FHUGE;
195	–	else
196	–	arad = (ndivs+ns)/arad;
197	–	if (arad > maxarad)
472		arad = maxarad;
473	<	else if (arad < minarad)
473	>	else
474	>	arad = (divcnt+hemi.ns)/arad;
475	>	if (pg != NULL) {               /* reduce radius if gradient large */
476	>	d = DOT(pg,pg);
477	>	if (d*arad*arad > 1.0)
478	>	arad = 1.0/sqrt(d);
479	>	}
480	>	if (arad < minarad) {
481		arad = minarad;
482	<	arad /= sqrt(r->rweight);
482	>	if (pg != NULL && d*arad*arad > 1.0) {  /* cap gradient */
483	>	d = 1.0/arad/sqrt(d);
484	>	for (i = 0; i < 3; i++)
485	>	pg[i] *= d;
486	>	}
487	>	}
488	>	if ((arad /= sqrt(wt)) > maxarad)
489	>	arad = maxarad;
490		return(arad);
203	–	oopsy:
204	–	if (div != NULL)
205	–	free((char *)div);
206	–	return(0.0);
491		}
492
493
494	<	inithemi(hp, r)                 /* initialize sampling hemisphere */
495	<	register AMBHEMI  *hp;
496	<	RAY  *r;
494	>	void
495	>	comperrs(                       /* compute initial error estimates */
496	>	AMBSAMP  *da,   /* assumes standard ordering */
497	>	AMBHEMI  *hp
498	>	)
499		{
214	–	register int  k;
215	–	/* set number of divisions */
216	–	hp->nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5;
217	–	hp->np = 2 * hp->nt;
218	–	/* make axes */
219	–	VCOPY(hp->uz, r->ron);
220	–	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
221	–	for (k = 0; k < 3; k++)
222	–	if (hp->uz[k] < 0.6 && hp->uz[k] > -0.6)
223	–	break;
224	–	if (k >= 3)
225	–	error(CONSISTENCY, "bad ray direction in inithemi");
226	–	hp->uy[k] = 1.0;
227	–	fcross(hp->ux, hp->uz, hp->uy);
228	–	normalize(hp->ux);
229	–	fcross(hp->uy, hp->ux, hp->uz);
230	–	}
231	–
232	–
233	–	comperrs(da, hp)                /* compute initial error estimates */
234	–	AMBSAMP  *da;
235	–	register AMBHEMI  *hp;
236	–	{
500		double  b, b2;
501		int  i, j;
502	<	register AMBSAMP  *dp;
502	>	AMBSAMP  *dp;
503		/* sum differences from neighbors */
504		dp = da;
505		for (i = 0; i < hp->nt; i++)
506		for (j = 0; j < hp->np; j++) {
507	+	#ifdef  DEBUG
508	+	if (dp->t != i || dp->p != j)
509	+	error(CONSISTENCY,
510	+	"division order in comperrs");
511	+	#endif
512		b = bright(dp[0].v);
513		if (i > 0) {            /* from above */
514		b2 = bright(dp[-hp->np].v) - b;
#	Line 253 | Line 521 | register AMBHEMI  *hp;
521		b2 *= b2 * 0.25;
522		dp[0].k += b2;
523		dp[-1].k += b2;
524	<	}
525	<	if (j == hp->np-1) {    /* around */
258	<	b2 = bright(dp[-(hp->np-1)].v) - b;
524	>	} else {                /* around */
525	>	b2 = bright(dp[hp->np-1].v) - b;
526		b2 *= b2 * 0.25;
527		dp[0].k += b2;
528	<	dp[-(hp->np-1)].k += b2;
528	>	dp[hp->np-1].k += b2;
529		}
530		dp++;
531		}
#	Line 276 | Line 543 | register AMBHEMI  *hp;
543		}
544
545
546	<	posgradient(gv, da, hp)                         /* compute position gradient */
547	<	FVECT  gv;
548	<	AMBSAMP  *da;
549	<	AMBHEMI  *hp;
546	>	void
547	>	posgradient(                                    /* compute position gradient */
548	>	FVECT  gv,
549	>	AMBSAMP  *da,                   /* assumes standard ordering */
550	>	AMBHEMI  *hp
551	>	)
552		{
553	<	gv[0] = 0.0; gv[1] = 0.0; gv[2] = 0.0;
553	>	int  i, j;
554	>	double  nextsine, lastsine, b, d;
555	>	double  mag0, mag1;
556	>	double  phi, cosp, sinp, xd, yd;
557	>	AMBSAMP  *dp;
558	>
559	>	xd = yd = 0.0;
560	>	for (j = 0; j < hp->np; j++) {
561	>	dp = da + j;
562	>	mag0 = mag1 = 0.0;
563	>	lastsine = 0.0;
564	>	for (i = 0; i < hp->nt; i++) {
565	>	#ifdef  DEBUG
566	>	if (dp->t != i || dp->p != j)
567	>	error(CONSISTENCY,
568	>	"division order in posgradient");
569	>	#endif
570	>	b = bright(dp->v);
571	>	if (i > 0) {
572	>	d = dp[-hp->np].r;
573	>	if (dp[0].r > d) d = dp[0].r;
574	>	/* sin(t)*cos(t)^2 */
575	>	d *= lastsine * (1.0 - (double)i/hp->nt);
576	>	mag0 += d*(b - bright(dp[-hp->np].v));
577	>	}
578	>	nextsine = sqrt((double)(i+1)/hp->nt);
579	>	if (j > 0) {
580	>	d = dp[-1].r;
581	>	if (dp[0].r > d) d = dp[0].r;
582	>	mag1 += d * (nextsine - lastsine) *
583	>	(b - bright(dp[-1].v));
584	>	} else {
585	>	d = dp[hp->np-1].r;
586	>	if (dp[0].r > d) d = dp[0].r;
587	>	mag1 += d * (nextsine - lastsine) *
588	>	(b - bright(dp[hp->np-1].v));
589	>	}
590	>	dp += hp->np;
591	>	lastsine = nextsine;
592	>	}
593	>	mag0 *= 2.0*PI / hp->np;
594	>	phi = 2.0*PI * (double)j/hp->np;
595	>	cosp = tcos(phi); sinp = tsin(phi);
596	>	xd += mag0*cosp - mag1*sinp;
597	>	yd += mag0*sinp + mag1*cosp;
598	>	}
599	>	for (i = 0; i < 3; i++)
600	>	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI;
601		}
602
603
604	<	dirgradient(gv, da, hp)                         /* compute direction gradient */
605	<	FVECT  gv;
606	<	AMBSAMP  *da;
607	<	AMBHEMI  *hp;
604	>	void
605	>	dirgradient(                                    /* compute direction gradient */
606	>	FVECT  gv,
607	>	AMBSAMP  *da,                   /* assumes standard ordering */
608	>	AMBHEMI  *hp
609	>	)
610		{
611	<	gv[0] = 0.0; gv[1] = 0.0; gv[2] = 0.0;
611	>	int  i, j;
612	>	double  mag;
613	>	double  phi, xd, yd;
614	>	AMBSAMP  *dp;
615	>
616	>	xd = yd = 0.0;
617	>	for (j = 0; j < hp->np; j++) {
618	>	dp = da + j;
619	>	mag = 0.0;
620	>	for (i = 0; i < hp->nt; i++) {
621	>	#ifdef  DEBUG
622	>	if (dp->t != i || dp->p != j)
623	>	error(CONSISTENCY,
624	>	"division order in dirgradient");
625	>	#endif
626	>	/* tan(t) */
627	>	mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
628	>	dp += hp->np;
629	>	}
630	>	phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
631	>	xd += mag * tcos(phi);
632	>	yd += mag * tsin(phi);
633	>	}
634	>	for (i = 0; i < 3; i++)
635	>	gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
636		}
637	+
638	+	#endif  /* ! NEWAMB */

Diff Legend

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