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.18 by greg, Mon Jun 6 19:14:28 2005 UTC vs.
Revision 2.26 by greg, Wed Apr 16 20:32:00 2014 UTC

#	Line 10 | Line 10 | static const char      RCSid[] = "$Id$";
10		#include "copyright.h"
11
12		#include  "ray.h"
13	–
13		#include  "ambient.h"
15	–
14		#include  "random.h"
15
16	+	#ifdef NEWAMB
17
18	+	extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
19	+
20	+	typedef struct {
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	+	#define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
32	+
33	+
34	+	static AMBHEMI *
35	+	inithemi(                       /* initialize sampling hemisphere */
36	+	COLOR   ac,
37	+	RAY     *r,
38	+	double  wt
39	+	)
40	+	{
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	+	ambsample(                              /* sample an ambient direction */
81	+	AMBHEMI *hp,
82	+	int     i,
83	+	int     j,
84	+	)
85	+	{
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	+
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	+	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	<	register AMBHEMI  *hp,
221	>	AMBHEMI  *hp,
222		COLOR ac,
223		RAY  *r,
224		double  wt
225		)
226		{
227		double  d;
228	<	register int  i;
228	>	int  i;
229		/* set number of divisions */
230		if (ambacc <= FTINY &&
231	<	wt > (d = 0.8*bright(ac)*r->rweight/(ambdiv*minweight)))
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 */
#	Line 58 | Line 258 | inithemi(                      /* initialize sampling hemisphere */
258
259		int
260		divsample(                              /* sample a division */
261	<	register AMBSAMP  *dp,
261	>	AMBSAMP  *dp,
262		AMBHEMI  *h,
263		RAY  *r
264		)
#	Line 69 | Line 269 | divsample(                             /* sample a division */
269		double  xd, yd, zd;
270		double  b2;
271		double  phi;
272	<	register int  i;
272	>	int  i;
273		/* ambient coefficient for weight */
274		if (ambacc > FTINY)
275		setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
#	Line 94 | Line 294 | divsample(                             /* sample a division */
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--;
#	Line 138 | Line 339 | ambnorm(                               /* standard order */
339		{
340		const AMBSAMP   *d1 = (const AMBSAMP *)p1;
341		const AMBSAMP   *d2 = (const AMBSAMP *)p2;
342	<	register int    c;
342	>	int     c;
343
344		if ( (c = d1->t - d2->t) )
345		return(c);
#	Line 148 | Line 349 | ambnorm(                               /* standard order */
349
350		double
351		doambient(                              /* compute ambient component */
352	<	COLOR  acol,
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;
367	<	register int  i, j;
366	>	int  divcnt;
367	>	int  i, j;
368		/* initialize hemisphere */
369	<	inithemi(&hemi, acol, r, wt);
370	<	ndivs = hemi.nt * hemi.np;
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(acol, 0.0, 0.0, 0.0);
377	<	if (ndivs == 0)
376	>	setcolor(rcol, 0.0, 0.0, 0.0);
377	>	if (divcnt == 0)
378		return(0.0);
379		/* allocate super-samples */
380		if (hemi.ns > 0 || pg != NULL || dg != NULL) {
381	<	div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
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;
#	Line 192 | Line 396 | doambient(                             /* compute ambient component */
396		dp->r = 0.0;
397		dp->n = 0;
398		if (divsample(dp, &hemi, r) < 0) {
399	<	if (div == NULL) continue;
400	<	dp++;
197	<	hemi.ns = 0;    /* incomplete sampling */
198	<	pg = dg = NULL;
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 (hemi.ns > 0 && arad > FTINY && ndivs/arad < minarad)
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 */
420	>	} else if (hemi.ns > 0) {       /* else perform super-sampling? */
421		comperrs(div, &hemi);                   /* compute errors */
422	<	qsort(div, ndivs, sizeof(AMBSAMP), ambcmp);     /* sort divs */
422	>	qsort(div, divcnt, sizeof(AMBSAMP), ambcmp);    /* sort divs */
423		/* super-sample */
424		for (i = hemi.ns; i > 0; i--) {
425		dnew = *div;
#	Line 217 | Line 428 | doambient(                             /* compute ambient component */
428		continue;
429		}
430		dp = div;               /* reinsert */
431	<	j = ndivs < i ? ndivs : i;
431	>	j = divcnt < i ? divcnt : i;
432		while (--j > 0 && dnew.k < dp[1].k) {
433		*dp = *(dp+1);
434		dp++;
#	Line 225 | Line 436 | doambient(                             /* compute ambient component */
436		*dp = dnew;
437		}
438		if (pg != NULL || dg != NULL)   /* restore order */
439	<	qsort(div, ndivs, sizeof(AMBSAMP), ambnorm);
439	>	qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
440		}
441		/* compute returned values */
442		if (div != NULL) {
443	<	arad = 0.0;
444	<	for (i = ndivs, dp = div; i-- > 0; dp++) {
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;
#	Line 256 | Line 467 | doambient(                             /* compute ambient component */
467		}
468		free((void *)div);
469		}
470	+	copycolor(rcol, acol);
471		if (arad <= FTINY)
472		arad = maxarad;
473		else
474	<	arad = (ndivs+hemi.ns)/arad;
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)
#	Line 282 | Line 494 | doambient(                             /* compute ambient component */
494		void
495		comperrs(                       /* compute initial error estimates */
496		AMBSAMP  *da,   /* assumes standard ordering */
497	<	register AMBHEMI  *hp
497	>	AMBHEMI  *hp
498		)
499		{
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++)
#	Line 335 | Line 547 | void
547		posgradient(                                    /* compute position gradient */
548		FVECT  gv,
549		AMBSAMP  *da,                   /* assumes standard ordering */
550	<	register AMBHEMI  *hp
550	>	AMBHEMI  *hp
551		)
552		{
553	<	register int  i, j;
553	>	int  i, j;
554		double  nextsine, lastsine, b, d;
555		double  mag0, mag1;
556		double  phi, cosp, sinp, xd, yd;
557	<	register AMBSAMP  *dp;
557	>	AMBSAMP  *dp;
558
559		xd = yd = 0.0;
560		for (j = 0; j < hp->np; j++) {
#	Line 393 | Line 605 | void
605		dirgradient(                                    /* compute direction gradient */
606		FVECT  gv,
607		AMBSAMP  *da,                   /* assumes standard ordering */
608	<	register AMBHEMI  *hp
608	>	AMBHEMI  *hp
609		)
610		{
611	<	register int  i, j;
611	>	int  i, j;
612		double  mag;
613		double  phi, xd, yd;
614	<	register AMBSAMP  *dp;
614	>	AMBSAMP  *dp;
615
616		xd = yd = 0.0;
617		for (j = 0; j < hp->np; j++) {
#	Line 422 | Line 634 | dirgradient(                                   /* compute direction gradient */
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)