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.58 by greg, Sun May 11 19:03:37 2014 UTC vs.
Revision 2.81 by greg, Thu Apr 12 18:02:45 2018 UTC

#	Line 21 | Line 21 | static const char      RCSid[] = "$Id$";
21		#include  "ambient.h"
22		#include  "random.h"
23
24	<	#ifdef NEWAMB
24	>	#ifndef OLDAMB
25
26		extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
27
#	Line 33 | Line 33 | typedef struct {
33
34		typedef struct {
35		RAY     *rp;            /* originating ray sample */
36	–	FVECT   ux, uy;         /* tangent axis unit vectors */
36		int     ns;             /* number of samples per axis */
37	+	int     sampOK;         /* acquired full sample set? */
38		COLOR   acoef;          /* division contribution coefficient */
39	+	double  acol[3];        /* accumulated color */
40	+	FVECT   ux, uy;         /* tangent axis unit vectors */
41		AMBSAMP sa[1];          /* sample array (extends struct) */
42		}  AMBHEMI;             /* ambient sample hemisphere */
43
#	Line 48 | Line 50 | typedef struct {
50		} FFTRI;                /* vectors and coefficients for Hessian calculation */
51
52
53	<	static AMBHEMI *
54	<	inithemi(                       /* initialize sampling hemisphere */
55	<	COLOR   ac,
56	<	RAY     *r,
57	<	double  wt
53	>	static int
54	>	ambcollision(                           /* proposed direciton collides? */
55	>	AMBHEMI *hp,
56	>	int     i,
57	>	int     j,
58	>	FVECT   dv
59		)
60		{
61	<	AMBHEMI *hp;
62	<	double  d;
63	<	int     n, i;
64	<	/* set number of divisions */
65	<	if (ambacc <= FTINY &&
66	<	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
67	<	wt = d;                 /* avoid ray termination */
68	<	n = sqrt(ambdiv * wt) + 0.5;
69	<	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
70	<	if (n < i)
71	<	n = i;
72	<	/* allocate sampling array */
73	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
74	<	if (hp == NULL)
75	<	return(NULL);
76	<	hp->rp = r;
77	<	hp->ns = n;
78	<	/* assign coefficient */
79	<	copycolor(hp->acoef, ac);
80	<	d = 1.0/(n*n);
81	<	scalecolor(hp->acoef, d);
82	<	/* make tangent plane axes */
83	<	hp->uy[0] = 0.5 - frandom();
84	<	hp->uy[1] = 0.5 - frandom();
85	<	hp->uy[2] = 0.5 - frandom();
86	<	for (i = 3; i--; )
84	<	if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
85	<	break;
86	<	if (i < 0)
87	<	error(CONSISTENCY, "bad ray direction in inithemi");
88	<	hp->uy[i] = 1.0;
89	<	VCROSS(hp->ux, hp->uy, r->ron);
90	<	normalize(hp->ux);
91	<	VCROSS(hp->uy, r->ron, hp->ux);
92	<	/* we're ready to sample */
93	<	return(hp);
61	>	double  cos_thresh;
62	>	int     ii, jj;
63	>	/* min. spacing = 1/4th division */
64	>	cos_thresh = (PI/4.)/(double)hp->ns;
65	>	cos_thresh = 1. - .5*cos_thresh*cos_thresh;
66	>	/* check existing neighbors */
67	>	for (ii = i-1; ii <= i+1; ii++) {
68	>	if (ii < 0) continue;
69	>	if (ii >= hp->ns) break;
70	>	for (jj = j-1; jj <= j+1; jj++) {
71	>	AMBSAMP *ap;
72	>	FVECT   avec;
73	>	double  dprod;
74	>	if (jj < 0) continue;
75	>	if (jj >= hp->ns) break;
76	>	if ((ii==i) & (jj==j)) continue;
77	>	ap = &ambsam(hp,ii,jj);
78	>	if (ap->d <= .5/FHUGE)
79	>	continue;       /* no one home */
80	>	VSUB(avec, ap->p, hp->rp->rop);
81	>	dprod = DOT(avec, dv);
82	>	if (dprod >= cos_thresh*VLEN(avec))
83	>	return(1);      /* collision */
84	>	}
85	>	}
86	>	return(0);                      /* nothing to worry about */
87		}
88
89
97	–	/* Sample ambient division and apply weighting coefficient */
90		static int
91	<	getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
91	>	ambsample(                              /* initial ambient division sample */
92	>	AMBHEMI *hp,
93	>	int     i,
94	>	int     j,
95	>	int     n
96	>	)
97		{
98	+	AMBSAMP *ap = &ambsam(hp,i,j);
99	+	RAY     ar;
100		int     hlist[3], ii;
101		double  spt[2], zd;
102	+	/* generate hemispherical sample */
103		/* ambient coefficient for weight */
104		if (ambacc > FTINY)
105	<	setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL);
105	>	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
106		else
107	<	copycolor(arp->rcoef, hp->acoef);
108	<	if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0)
107	>	copycolor(ar.rcoef, hp->acoef);
108	>	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
109		return(0);
110		if (ambacc > FTINY) {
111	<	multcolor(arp->rcoef, hp->acoef);
112	<	scalecolor(arp->rcoef, 1./AVGREFL);
111	>	multcolor(ar.rcoef, hp->acoef);
112	>	scalecolor(ar.rcoef, 1./AVGREFL);
113		}
114		hlist[0] = hp->rp->rno;
115		hlist[1] = j;
116		hlist[2] = i;
117		multisamp(spt, 2, urand(ilhash(hlist,3)+n));
118	<	if (!n) {                       /* avoid border samples for n==0 */
119	<	if ((spt[0] < 0.1) | (spt[0] >= 0.9))
120	<	spt[0] = 0.1 + 0.8*frandom();
121	<	if ((spt[1] < 0.1) | (spt[1] >= 0.9))
122	<	spt[1] = 0.1 + 0.8*frandom();
123	<	}
118	>	resample:
119		SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
120		zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
121		for (ii = 3; ii--; )
122	<	arp->rdir[ii] = spt[0]*hp->ux[ii] +
122	>	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
123		spt[1]*hp->uy[ii] +
124		zd*hp->rp->ron[ii];
125	<	checknorm(arp->rdir);
125	>	checknorm(ar.rdir);
126	>	/* avoid coincident samples */
127	>	if (!n && ambcollision(hp, i, j, ar.rdir)) {
128	>	spt[0] = frandom(); spt[1] = frandom();
129	>	goto resample;          /* reject this sample */
130	>	}
131		dimlist[ndims++] = AI(hp,i,j) + 90171;
132	<	rayvalue(arp);                  /* evaluate ray */
133	<	ndims--;                        /* apply coefficient */
134	<	multcolor(arp->rcol, arp->rcoef);
132	>	rayvalue(&ar);                  /* evaluate ray */
133	>	ndims--;
134	>	if (ar.rt <= FTINY)
135	>	return(0);              /* should never happen */
136	>	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
137	>	if (ar.rt*ap->d < 1.0)          /* new/closer distance? */
138	>	ap->d = 1.0/ar.rt;
139	>	if (!n) {                       /* record first vertex & value */
140	>	if (ar.rt > 10.0*thescene.cusize + 1000.)
141	>	ar.rt = 10.0*thescene.cusize + 1000.;
142	>	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
143	>	copycolor(ap->v, ar.rcol);
144	>	} else {                        /* else update recorded value */
145	>	hp->acol[RED] -= colval(ap->v,RED);
146	>	hp->acol[GRN] -= colval(ap->v,GRN);
147	>	hp->acol[BLU] -= colval(ap->v,BLU);
148	>	zd = 1.0/(double)(n+1);
149	>	scalecolor(ar.rcol, zd);
150	>	zd *= (double)n;
151	>	scalecolor(ap->v, zd);
152	>	addcolor(ap->v, ar.rcol);
153	>	}
154	>	addcolor(hp->acol, ap->v);      /* add to our sum */
155		return(1);
156		}
157
158
159	<	static AMBSAMP *
140	<	ambsample(                              /* initial ambient division sample */
141	<	AMBHEMI *hp,
142	<	int     i,
143	<	int     j
144	<	)
145	<	{
146	<	AMBSAMP *ap = &ambsam(hp,i,j);
147	<	RAY     ar;
148	<	/* generate hemispherical sample */
149	<	if (!getambsamp(&ar, hp, i, j, 0) || ar.rt <= FTINY) {
150	<	memset(ap, 0, sizeof(AMBSAMP));
151	<	return(NULL);
152	<	}
153	<	ap->d = 1.0/ar.rt;              /* limit vertex distance */
154	<	if (ar.rt > 10.0*thescene.cusize)
155	<	ar.rt = 10.0*thescene.cusize;
156	<	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
157	<	copycolor(ap->v, ar.rcol);
158	<	return(ap);
159	<	}
160	<
161	<
162	<	/* Estimate errors based on ambient division differences */
159	>	/* Estimate variance based on relative ambient division differences */
160		static float *
161		getambdiffs(AMBHEMI *hp)
162		{
163	+	const double    normf = 1./bright(hp->acoef);
164		float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
165		float   *ep;
166		AMBSAMP *ap;
167	<	double  b, d2;
167	>	double  b, b1, d2;
168		int     i, j;
169
170		if (earr == NULL)               /* out of memory? */
171		return(NULL);
172	<	/* compute squared neighbor diffs */
172	>	/* sum squared neighbor diffs */
173		for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
174		for (j = 0; j < hp->ns; j++, ap++, ep++) {
175		b = bright(ap[0].v);
176		if (i) {                /* from above */
177	<	d2 = b - bright(ap[-hp->ns].v);
178	<	d2 *= d2;
177	>	b1 = bright(ap[-hp->ns].v);
178	>	d2 = (b - b1)/(b + b1);
179	>	d2 *= d2*normf;
180		ep[0] += d2;
181		ep[-hp->ns] += d2;
182		}
183		if (!j) continue;
184		/* from behind */
185	<	d2 = b - bright(ap[-1].v);
186	<	d2 *= d2;
185	>	b1 = bright(ap[-1].v);
186	>	d2 = (b - b1)/(b + b1);
187	>	d2 *= d2*normf;
188		ep[0] += d2;
189		ep[-1] += d2;
190		if (!i) continue;
191		/* diagonal */
192	<	d2 = b - bright(ap[-hp->ns-1].v);
193	<	d2 *= d2;
192	>	b1 = bright(ap[-hp->ns-1].v);
193	>	d2 = (b - b1)/(b + b1);
194	>	d2 *= d2*normf;
195		ep[0] += d2;
196		ep[-hp->ns-1] += d2;
197		}
#	Line 213 | Line 214 | getambdiffs(AMBHEMI *hp)
214
215		/* Perform super-sampling on hemisphere (introduces bias) */
216		static void
217	<	ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
217	>	ambsupersamp(AMBHEMI *hp, int cnt)
218		{
219		float   *earr = getambdiffs(hp);
220		double  e2rem = 0;
220	–	AMBSAMP *ap;
221	–	RAY     ar;
222	–	double  asum[3];
221		float   *ep;
222		int     i, j, n, nss;
223
#	Line 229 | Line 227 | ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
227		for (ep = earr + hp->ns*hp->ns; ep > earr; )
228		e2rem += *--ep;
229		ep = earr;                      /* perform super-sampling */
230	<	for (ap = hp->sa, i = 0; i < hp->ns; i++)
231	<	for (j = 0; j < hp->ns; j++, ap++) {
230	>	for (i = 0; i < hp->ns; i++)
231	>	for (j = 0; j < hp->ns; j++) {
232		if (e2rem <= FTINY)
233		goto done;      /* nothing left to do */
234		nss = *ep/e2rem*cnt + frandom();
235	<	asum[0] = asum[1] = asum[2] = 0.0;
236	<	for (n = 1; n <= nss; n++) {
237	<	if (!getambsamp(&ar, hp, i, j, n)) {
240	<	nss = n-1;
241	<	break;
242	<	}
243	<	addcolor(asum, ar.rcol);
244	<	}
245	<	if (nss) {              /* update returned ambient value */
246	<	const double    ssf = 1./(nss + 1.);
247	<	for (n = 3; n--; )
248	<	acol[n] += ssf*asum[n] +
249	<	(ssf - 1.)*colval(ap->v,n);
250	<	}
251	<	e2rem -= *ep++;         /* update remainders */
252	<	cnt -= nss;
235	>	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
236	>	if (!--cnt) goto done;
237	>	e2rem -= *ep++;         /* update remainder */
238		}
239		done:
240		free(earr);
241		}
242
243
244	+	static AMBHEMI *
245	+	samp_hemi(                              /* sample indirect hemisphere */
246	+	COLOR   rcol,
247	+	RAY     *r,
248	+	double  wt
249	+	)
250	+	{
251	+	AMBHEMI *hp;
252	+	double  d;
253	+	int     n, i, j;
254	+	/* insignificance check */
255	+	if (bright(rcol) <= FTINY)
256	+	return(NULL);
257	+	/* set number of divisions */
258	+	if (ambacc <= FTINY &&
259	+	wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight)))
260	+	wt = d;                 /* avoid ray termination */
261	+	n = sqrt(ambdiv * wt) + 0.5;
262	+	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
263	+	if (n < i)
264	+	n = i;
265	+	/* allocate sampling array */
266	+	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
267	+	if (hp == NULL)
268	+	error(SYSTEM, "out of memory in samp_hemi");
269	+	hp->rp = r;
270	+	hp->ns = n;
271	+	hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0;
272	+	memset(hp->sa, 0, sizeof(AMBSAMP)*n*n);
273	+	hp->sampOK = 0;
274	+	/* assign coefficient */
275	+	copycolor(hp->acoef, rcol);
276	+	d = 1.0/(n*n);
277	+	scalecolor(hp->acoef, d);
278	+	/* make tangent plane axes */
279	+	if (!getperpendicular(hp->ux, r->ron, 1))
280	+	error(CONSISTENCY, "bad ray direction in samp_hemi");
281	+	VCROSS(hp->uy, r->ron, hp->ux);
282	+	/* sample divisions */
283	+	for (i = hp->ns; i--; )
284	+	for (j = hp->ns; j--; )
285	+	hp->sampOK += ambsample(hp, i, j, 0);
286	+	copycolor(rcol, hp->acol);
287	+	if (!hp->sampOK) {              /* utter failure? */
288	+	free(hp);
289	+	return(NULL);
290	+	}
291	+	if (hp->sampOK < hp->ns*hp->ns) {
292	+	hp->sampOK *= -1;       /* soft failure */
293	+	return(hp);
294	+	}
295	+	n = ambssamp*wt + 0.5;
296	+	if (n > 8) {                    /* perform super-sampling? */
297	+	ambsupersamp(hp, n);
298	+	copycolor(rcol, hp->acol);
299	+	}
300	+	return(hp);                     /* all is well */
301	+	}
302	+
303	+
304		/* Return brightness of farthest ambient sample */
305		static double
306		back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
#	Line 596 | Line 641 | static uint32
641		ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
642		{
643		const double    max_d = 1.0/(minarad*ambacc + 0.001);
644	<	const double    ang_res = 0.5*PI/(hp->ns-1);
645	<	const double    ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY));
644	>	const double    ang_res = 0.5*PI/hp->ns;
645	>	const double    ang_step = ang_res/((int)(16/PI*ang_res) + 1.01);
646		double          avg_d = 0;
647		uint32          flgs = 0;
648		FVECT           vec;
#	Line 605 | Line 650 | ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, c
650		double          ang, a1;
651		int             i, j;
652		/* don't bother for a few samples */
653	<	if (hp->ns < 12)
653	>	if (hp->ns < 8)
654		return(0);
655		/* check distances overhead */
656		for (i = hp->ns*3/4; i-- > hp->ns>>2; )
#	Line 623 | Line 668 | ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, c
668		if ((ap->d <= FTINY) | (ap->d >= max_d))
669		continue;       /* too far or too near */
670		VSUB(vec, ap->p, hp->rp->rop);
671	<	u = DOT(vec, uv[0]) * ap->d;
672	<	v = DOT(vec, uv[1]) * ap->d;
673	<	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0)
671	>	u = DOT(vec, uv[0]);
672	>	v = DOT(vec, uv[1]);
673	>	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v)
674		continue;       /* occluder outside ellipse */
675		ang = atan2a(v, u);     /* else set direction flags */
676	<	for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step)
676	>	for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step)
677		flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0)));
678		}
634	–	/* add low-angle incident (< 20deg) */
635	–	if (fabs(hp->rp->rod) <= 0.342) {
636	–	u = -DOT(hp->rp->rdir, uv[0]);
637	–	v = -DOT(hp->rp->rdir, uv[1]);
638	–	if ((r0*r0*u*u + r1*r1*v*v) > hp->rp->rot*hp->rp->rot) {
639	–	ang = atan2a(v, u);
640	–	ang += 2.*PI*(ang < 0);
641	–	ang *= 16/PI;
642	–	if ((ang < .5) | (ang >= 31.5))
643	–	flgs |= 0x80000001;
644	–	else
645	–	flgs |= 3L<<(int)(ang-.5);
646	–	}
647	–	}
679		return(flgs);
680		}
681
#	Line 661 | Line 692 | doambient(                             /* compute ambient component */
692		uint32  *crlp                   /* returned (optional) */
693		)
694		{
695	<	AMBHEMI *hp = inithemi(rcol, r, wt);
665	<	int     cnt;
695	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
696		FVECT   my_uv[2];
697	<	double  d, K, acol[3];
697	>	double  d, K;
698		AMBSAMP *ap;
699	<	int     i, j;
700	<	/* check/initialize */
671	<	if (hp == NULL)
672	<	return(0);
699	>	int     i;
700	>	/* clear return values */
701		if (uv != NULL)
702		memset(uv, 0, sizeof(FVECT)*2);
703		if (ra != NULL)
#	Line 680 | Line 708 | doambient(                             /* compute ambient component */
708		dg[0] = dg[1] = 0.0;
709		if (crlp != NULL)
710		*crlp = 0;
711	<	/* sample the hemisphere */
712	<	acol[0] = acol[1] = acol[2] = 0.0;
713	<	cnt = 0;
714	<	for (i = hp->ns; i--; )
715	<	for (j = hp->ns; j--; )
716	<	if ((ap = ambsample(hp, i, j)) != NULL) {
717	<	addcolor(acol, ap->v);
690	<	++cnt;
691	<	}
692	<	if (!cnt) {
693	<	setcolor(rcol, 0.0, 0.0, 0.0);
694	<	free(hp);
695	<	return(0);              /* no valid samples */
711	>	if (hp == NULL)                 /* sampling falure? */
712	>	return(0);
713	>
714	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) ||
715	>	(hp->sampOK < 0) | (hp->ns < 6)) {
716	>	free(hp);               /* Hessian not requested/possible */
717	>	return(-1);             /* value-only return value */
718		}
719	<	if (cnt < hp->ns*hp->ns) {      /* incomplete sampling? */
698	<	copycolor(rcol, acol);
699	<	free(hp);
700	<	return(-1);             /* return value w/o Hessian */
701	<	}
702	<	cnt = ambssamp*wt + 0.5;        /* perform super-sampling? */
703	<	if (cnt > 8)
704	<	ambsupersamp(acol, hp, cnt);
705	<	copycolor(rcol, acol);          /* final indirect irradiance/PI */
706	<	if ((ra == NULL) & (pg == NULL) & (dg == NULL)) {
707	<	free(hp);
708	<	return(-1);             /* no radius or gradient calc. */
709	<	}
710	<	if ((d = bright(acol)) > FTINY) {       /* normalize Y values */
719	>	if ((d = bright(rcol)) > FTINY) {       /* normalize Y values */
720		d = 0.99*(hp->ns*hp->ns)/d;
721		K = 0.01;
722		} else {                        /* or fall back on geometric Hessian */
#	Line 742 | Line 751 | doambient(                             /* compute ambient component */
751		if (ra[1] < minarad)
752		ra[1] = minarad;
753		}
754	<	ra[0] *= d = 1.0/sqrt(sqrt(wt));
754	>	ra[0] *= d = 1.0/sqrt(wt);
755		if ((ra[1] *= d) > 2.0*ra[0])
756		ra[1] = 2.0*ra[0];
757		if (ra[1] > maxarad) {
#	Line 751 | Line 760 | doambient(                             /* compute ambient component */
760		ra[0] = maxarad;
761		}
762		/* flag encroached directions */
763	<	if ((wt >= 0.89*AVGREFL) & (crlp != NULL))
763	>	if (crlp != NULL)
764		*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc);
765		if (pg != NULL) {       /* cap gradient if necessary */
766		d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];

Diff Legend

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