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.11 by greg, Tue Aug 13 12:16:41 1991 UTC vs.
Revision 2.9 by greg, Sat Feb 22 02:07:28 2003 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	+	/* ====================================================================
11	+	* The Radiance Software License, Version 1.0
12	+	*
13	+	* Copyright (c) 1990 - 2002 The Regents of the University of California,
14	+	* through Lawrence Berkeley National Laboratory.   All rights reserved.
15	+	*
16	+	* Redistribution and use in source and binary forms, with or without
17	+	* modification, are permitted provided that the following conditions
18	+	* are met:
19	+	*
20	+	* 1. Redistributions of source code must retain the above copyright
21	+	*         notice, this list of conditions and the following disclaimer.
22	+	*
23	+	* 2. Redistributions in binary form must reproduce the above copyright
24	+	*       notice, this list of conditions and the following disclaimer in
25	+	*       the documentation and/or other materials provided with the
26	+	*       distribution.
27	+	*
28	+	* 3. The end-user documentation included with the redistribution,
29	+	*           if any, must include the following acknowledgment:
30	+	*             "This product includes Radiance software
31	+	*                 (http://radsite.lbl.gov/)
32	+	*                 developed by the Lawrence Berkeley National Laboratory
33	+	*               (http://www.lbl.gov/)."
34	+	*       Alternately, this acknowledgment may appear in the software itself,
35	+	*       if and wherever such third-party acknowledgments normally appear.
36	+	*
37	+	* 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
38	+	*       and "The Regents of the University of California" must
39	+	*       not be used to endorse or promote products derived from this
40	+	*       software without prior written permission. For written
41	+	*       permission, please contact [email protected].
42	+	*
43	+	* 5. Products derived from this software may not be called "Radiance",
44	+	*       nor may "Radiance" appear in their name, without prior written
45	+	*       permission of Lawrence Berkeley National Laboratory.
46	+	*
47	+	* THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
48	+	* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
49	+	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
50	+	* DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
51	+	* ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
52	+	* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
53	+	* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
54	+	* USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
55	+	* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
56	+	* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
57	+	* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
58	+	* SUCH DAMAGE.
59	+	* ====================================================================
60	+	*
61	+	* This software consists of voluntary contributions made by many
62	+	* individuals on behalf of Lawrence Berkeley National Laboratory.   For more
63	+	* information on Lawrence Berkeley National Laboratory, please see
64	+	* <http://www.lbl.gov/>.
65	+	*/
66	+
67		#include  "ray.h"
68
69		#include  "ambient.h"
70
71		#include  "random.h"
72
17	–	typedef struct {
18	–	short  t, p;            /* theta, phi indices */
19	–	COLOR  v;               /* value sum */
20	–	float  r;               /* 1/distance sum */
21	–	float  k;               /* variance for this division */
22	–	int  n;                 /* number of subsamples */
23	–	}  AMBSAMP;             /* ambient sample division */
73
25	–	typedef struct {
26	–	FVECT  ux, uy, uz;      /* x, y and z axis directions */
27	–	short  nt, np;          /* number of theta and phi directions */
28	–	}  AMBHEMI;             /* ambient sample hemisphere */
29	–
30	–	extern double  sin(), cos(), sqrt();
31	–
32	–
74		static int
75		ambcmp(d1, d2)                          /* decreasing order */
76		AMBSAMP  *d1, *d2;
#	Line 54 | Line 95 | AMBSAMP  *d1, *d2;
95		}
96
97
98	+	int
99		divsample(dp, h, r)                     /* sample a division */
100		register AMBSAMP  *dp;
101		AMBHEMI  *h;
#	Line 67 | Line 109 | RAY  *r;
109		double  phi;
110		register int  i;
111
112	<	if (rayorigin(&ar, r, AMBIENT, 0.5) < 0)
112	>	if (rayorigin(&ar, r, AMBIENT, AVGREFL) < 0)
113		return(-1);
114		hlist[0] = r->rno;
115		hlist[1] = dp->t;
116		hlist[2] = dp->p;
117	<	peano(spt, 2, urand(ilhash(hlist,3)+dp->n), .01);
117	>	multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
118		zd = sqrt((dp->t + spt[0])/h->nt);
119		phi = 2.0*PI * (dp->p + spt[1])/h->np;
120	<	xd = cos(phi) * zd;
121	<	yd = sin(phi) * zd;
120	>	xd = tcos(phi) * zd;
121	>	yd = tsin(phi) * zd;
122		zd = sqrt(1.0 - zd*zd);
123		for (i = 0; i < 3; i++)
124		ar.rdir[i] =    xd*h->ux[i] +
#	Line 86 | Line 128 | RAY  *r;
128		rayvalue(&ar);
129		ndims--;
130		addcolor(dp->v, ar.rcol);
131	+	/* use rt to improve gradient calc */
132		if (ar.rt > FTINY && ar.rt < FHUGE)
133		dp->r += 1.0/ar.rt;
134		/* (re)initialize error */
#	Line 100 | Line 143 | RAY  *r;
143
144
145		double
146	<	doambient(acol, r, pg, dg)              /* compute ambient component */
146	>	doambient(acol, r, wt, pg, dg)          /* compute ambient component */
147		COLOR  acol;
148		RAY  *r;
149	+	double  wt;
150		FVECT  pg, dg;
151		{
152		double  b, d;
#	Line 116 | Line 160 | FVECT  pg, dg;
160		/* initialize color */
161		setcolor(acol, 0.0, 0.0, 0.0);
162		/* initialize hemisphere */
163	<	inithemi(&hemi, r);
163	>	inithemi(&hemi, r, wt);
164		ndivs = hemi.nt * hemi.np;
165		if (ndivs == 0)
166		return(0.0);
167		/* set number of super-samples */
168	<	ns = ambssamp * r->rweight + 0.5;
168	>	ns = ambssamp * wt + 0.5;
169		if (ns > 0 || pg != NULL || dg != NULL) {
170		div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
171		if (div == NULL)
#	Line 140 | Line 184 | FVECT  pg, dg;
184		dp->n = 0;
185		if (divsample(dp, &hemi, r) < 0)
186		goto oopsy;
187	+	arad += dp->r;
188		if (div != NULL)
189		dp++;
190	<	else {
190	>	else
191		addcolor(acol, dp->v);
147	–	arad += dp->r;
148	–	}
192		}
193	<	if (ns > 0) {                   /* perform super-sampling */
193	>	if (ns > 0 && arad > FTINY && ndivs/arad < minarad)
194	>	ns = 0;                 /* close enough */
195	>	else if (ns > 0) {              /* else perform super-sampling */
196		comperrs(div, &hemi);                   /* compute errors */
197		qsort(div, ndivs, sizeof(AMBSAMP), ambcmp);     /* sort divs */
198		/* super-sample */
#	Line 169 | Line 214 | FVECT  pg, dg;
214		}
215		/* compute returned values */
216		if (div != NULL) {
217	+	arad = 0.0;
218		for (i = ndivs, dp = div; i-- > 0; dp++) {
219		arad += dp->r;
220		if (dp->n > 1) {
#	Line 200 | Line 246 | FVECT  pg, dg;
246		for (i = 0; i < 3; i++)
247		dg[i] = 0.0;
248		}
249	<	free((char *)div);
249	>	free((void *)div);
250		}
251		b = 1.0/ndivs;
252		scalecolor(acol, b);
253		if (arad <= FTINY)
254	<	arad = FHUGE;
254	>	arad = maxarad;
255		else
256		arad = (ndivs+ns)/arad;
257	<	if (arad > maxarad)
258	<	arad = maxarad;
259	<	else if (arad < minarad)
257	>	if (pg != NULL) {               /* reduce radius if gradient large */
258	>	d = DOT(pg,pg);
259	>	if (d*arad*arad > 1.0)
260	>	arad = 1.0/sqrt(d);
261	>	}
262	>	if (arad < minarad) {
263		arad = minarad;
264	<	arad /= sqrt(r->rweight);
265	<	if (pg != NULL) {               /* clip pos. gradient if too large */
217	<	d = 4.0*DOT(pg,pg)*arad*arad;
218	<	if (d > 1.0) {
219	<	d = 1.0/sqrt(d);
264	>	if (pg != NULL && d*arad*arad > 1.0) {  /* cap gradient */
265	>	d = 1.0/arad/sqrt(d);
266		for (i = 0; i < 3; i++)
267		pg[i] *= d;
268		}
269		}
270	+	if ((arad /= sqrt(wt)) > maxarad)
271	+	arad = maxarad;
272		return(arad);
273		oopsy:
274		if (div != NULL)
275	<	free((char *)div);
275	>	free((void *)div);
276		return(0.0);
277		}
278
279
280	<	inithemi(hp, r)                 /* initialize sampling hemisphere */
280	>	void
281	>	inithemi(hp, r, wt)             /* initialize sampling hemisphere */
282		register AMBHEMI  *hp;
283		RAY  *r;
284	+	double  wt;
285		{
286		register int  i;
287		/* set number of divisions */
288	<	hp->nt = sqrt(ambdiv * r->rweight / PI) + 0.5;
289	<	hp->np = PI * hp->nt;
288	>	if (wt < (.25*PI)/ambdiv+FTINY) {
289	>	hp->nt = hp->np = 0;
290	>	return;                 /* zero samples */
291	>	}
292	>	hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
293	>	hp->np = PI * hp->nt + 0.5;
294		/* make axes */
295		VCOPY(hp->uz, r->ron);
296		hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
#	Line 252 | Line 306 | RAY  *r;
306		}
307
308
309	+	void
310		comperrs(da, hp)                /* compute initial error estimates */
311		AMBSAMP  *da;           /* assumes standard ordering */
312		register AMBHEMI  *hp;
#	Line 302 | Line 357 | register AMBHEMI  *hp;
357		}
358
359
360	+	void
361		posgradient(gv, da, hp)                         /* compute position gradient */
362		FVECT  gv;
363		AMBSAMP  *da;                   /* assumes standard ordering */
364	<	AMBHEMI  *hp;
364	>	register AMBHEMI  *hp;
365		{
366		register int  i, j;
367	<	double  b, d;
367	>	double  nextsine, lastsine, b, d;
368		double  mag0, mag1;
369		double  phi, cosp, sinp, xd, yd;
370		register AMBSAMP  *dp;
#	Line 317 | Line 373 | AMBHEMI  *hp;
373		for (j = 0; j < hp->np; j++) {
374		dp = da + j;
375		mag0 = mag1 = 0.0;
376	+	lastsine = 0.0;
377		for (i = 0; i < hp->nt; i++) {
378		#ifdef  DEBUG
379		if (dp->t != i || dp->p != j)
#	Line 327 | Line 384 | AMBHEMI  *hp;
384		if (i > 0) {
385		d = dp[-hp->np].r;
386		if (dp[0].r > d) d = dp[0].r;
387	<	d *= 1.0 - sqrt((double)i/hp->nt);
387	>	/* sin(t)*cos(t)^2 */
388	>	d *= lastsine * (1.0 - (double)i/hp->nt);
389		mag0 += d*(b - bright(dp[-hp->np].v));
390		}
391	+	nextsine = sqrt((double)(i+1)/hp->nt);
392		if (j > 0) {
393		d = dp[-1].r;
394		if (dp[0].r > d) d = dp[0].r;
395	<	mag1 += d*(b - bright(dp[-1].v));
395	>	mag1 += d * (nextsine - lastsine) *
396	>	(b - bright(dp[-1].v));
397		} else {
398		d = dp[hp->np-1].r;
399		if (dp[0].r > d) d = dp[0].r;
400	<	mag1 += d*(b - bright(dp[hp->np-1].v));
400	>	mag1 += d * (nextsine - lastsine) *
401	>	(b - bright(dp[hp->np-1].v));
402		}
403		dp += hp->np;
404	+	lastsine = nextsine;
405		}
406	<	if (hp->nt > 1) {
345	<	mag0 /= (double)hp->np;
346	<	mag1 /= (double)hp->nt;
347	<	}
406	>	mag0 *= 2.0*PI / hp->np;
407		phi = 2.0*PI * (double)j/hp->np;
408	<	cosp = cos(phi); sinp = sin(phi);
408	>	cosp = tcos(phi); sinp = tsin(phi);
409		xd += mag0*cosp - mag1*sinp;
410		yd += mag0*sinp + mag1*cosp;
411		}
#	Line 355 | Line 414 | AMBHEMI  *hp;
414		}
415
416
417	+	void
418		dirgradient(gv, da, hp)                         /* compute direction gradient */
419		FVECT  gv;
420		AMBSAMP  *da;                   /* assumes standard ordering */
421	<	AMBHEMI  *hp;
421	>	register AMBHEMI  *hp;
422		{
423		register int  i, j;
424		double  mag;
#	Line 375 | Line 435 | AMBHEMI  *hp;
435		error(CONSISTENCY,
436		"division order in dirgradient");
437		#endif
438	<	mag += sqrt((i+.5)/hp->nt)*bright(dp->v);
438	>	/* tan(t) */
439	>	mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
440		dp += hp->np;
441		}
442		phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
443	<	xd += mag * cos(phi);
444	<	yd += mag * sin(phi);
443	>	xd += mag * tcos(phi);
444	>	yd += mag * tsin(phi);
445		}
446		for (i = 0; i < 3; i++)
447	<	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*PI/(hp->nt*hp->np);
447	>	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])/(hp->nt*hp->np);
448		}

Diff Legend

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