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.51 by greg, Wed May 7 20:20:24 2014 UTC vs.
Revision 2.79 by greg, Tue Jan 9 00:52:35 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
28	–	/* vertex direction bit positions */
29	–	#define VDB_xy  0
30	–	#define VDB_y   01
31	–	#define VDB_x   02
32	–	#define VDB_Xy  03
33	–	#define VDB_xY  04
34	–	#define VDB_X   05
35	–	#define VDB_Y   06
36	–	#define VDB_XY  07
37	–	/* get opposite vertex direction bit */
38	–	#define VDB_OPP(f)      (~(f) & 07)
39	–	/* adjacent triangle vertex flags */
40	–	static const int  adjacent_trifl[8] = {
41	–	0,                      /* forbidden diagonal */
42	–	1<<VDB_x|1<<VDB_y|1<<VDB_Xy,
43	–	1<<VDB_y|1<<VDB_x|1<<VDB_xY,
44	–	1<<VDB_y|1<<VDB_Xy|1<<VDB_X,
45	–	1<<VDB_x|1<<VDB_xY|1<<VDB_Y,
46	–	1<<VDB_Xy|1<<VDB_X|1<<VDB_Y,
47	–	1<<VDB_xY|1<<VDB_Y|1<<VDB_X,
48	–	0,                      /* forbidden diagonal */
49	–	};
50	–
28		typedef struct {
29		COLOR   v;              /* hemisphere sample value */
30		float   d;              /* reciprocal distance (1/rt) */
#	Line 56 | Line 33 | typedef struct {
33
34		typedef struct {
35		RAY     *rp;            /* originating ray sample */
59	–	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
44	<	#define ambndx(h,i,j)   ((i)*(h)->ns + (j))
45	<	#define ambsam(h,i,j)   (h)->sa[ambndx(h,i,j)]
44	>	#define AI(h,i,j)       ((i)*(h)->ns + (j))
45	>	#define ambsam(h,i,j)   (h)->sa[AI(h,i,j)]
46
47		typedef struct {
48		FVECT   r_i, r_i1, e_i, rcp, rI2_eJ2;
49		double  I1, I2;
71	–	int     valid;
50		} FFTRI;                /* vectors and coefficients for Hessian calculation */
51
52
75	–	/* Get index for adjacent vertex */
53		static int
54	<	adjacent_verti(AMBHEMI *hp, int i, int j, int dbit)
54	>	ambcollision(                           /* proposed direciton collides? */
55	>	AMBHEMI *hp,
56	>	int     i,
57	>	int     j,
58	>	FVECT   dv
59	>	)
60		{
61	<	int     i0 = i*hp->ns + j;
62	<
63	<	switch (dbit) {
64	<	case VDB_y:     return(i0 - hp->ns);
65	<	case VDB_x:     return(i0 - 1);
66	<	case VDB_Xy:    return(i0 - hp->ns + 1);
67	<	case VDB_xY:    return(i0 + hp->ns - 1);
68	<	case VDB_X:     return(i0 + 1);
69	<	case VDB_Y:     return(i0 + hp->ns);
70	<	/* the following should never occur */
71	<	case VDB_xy:    return(i0 - hp->ns - 1);
72	<	case VDB_XY:    return(i0 + hp->ns + 1);
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(-1);
86	>	return(0);                      /* nothing to worry about */
87		}
88
89
96	–	/* Get vertex direction bit for the opposite edge to complete triangle */
90		static int
91	<	vdb_edge(int db1, int db2)
92	<	{
93	<	switch (db1) {
94	<	case VDB_x:     return(db2==VDB_y ? VDB_Xy : VDB_Y);
95	<	case VDB_y:     return(db2==VDB_x ? VDB_xY : VDB_X);
103	<	case VDB_X:     return(db2==VDB_Xy ? VDB_y : VDB_xY);
104	<	case VDB_Y:     return(db2==VDB_xY ? VDB_x : VDB_Xy);
105	<	case VDB_xY:    return(db2==VDB_x ? VDB_y : VDB_X);
106	<	case VDB_Xy:    return(db2==VDB_y ? VDB_x : VDB_Y);
107	<	}
108	<	error(INTERNAL, "forbidden diagonal in vdb_edge()");
109	<	return(-1);
110	<	}
111	<
112	<
113	<	static AMBHEMI *
114	<	inithemi(                       /* initialize sampling hemisphere */
115	<	COLOR   ac,
116	<	RAY     *r,
117	<	double  wt
91	>	ambsample(                              /* initial ambient division sample */
92	>	AMBHEMI *hp,
93	>	int     i,
94	>	int     j,
95	>	int     n
96		)
97		{
98	<	AMBHEMI *hp;
99	<	double  d;
122	<	int     n, i;
123	<	/* set number of divisions */
124	<	if (ambacc <= FTINY &&
125	<	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
126	<	wt = d;                 /* avoid ray termination */
127	<	n = sqrt(ambdiv * wt) + 0.5;
128	<	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
129	<	if (n < i)
130	<	n = i;
131	<	/* allocate sampling array */
132	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
133	<	if (hp == NULL)
134	<	return(NULL);
135	<	hp->rp = r;
136	<	hp->ns = n;
137	<	/* assign coefficient */
138	<	copycolor(hp->acoef, ac);
139	<	d = 1.0/(n*n);
140	<	scalecolor(hp->acoef, d);
141	<	/* make tangent plane axes */
142	<	hp->uy[0] = 0.5 - frandom();
143	<	hp->uy[1] = 0.5 - frandom();
144	<	hp->uy[2] = 0.5 - frandom();
145	<	for (i = 3; i--; )
146	<	if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
147	<	break;
148	<	if (i < 0)
149	<	error(CONSISTENCY, "bad ray direction in inithemi");
150	<	hp->uy[i] = 1.0;
151	<	VCROSS(hp->ux, hp->uy, r->ron);
152	<	normalize(hp->ux);
153	<	VCROSS(hp->uy, r->ron, hp->ux);
154	<	/* we're ready to sample */
155	<	return(hp);
156	<	}
157	<
158	<
159	<	/* Sample ambient division and apply weighting coefficient */
160	<	static int
161	<	getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
162	<	{
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 */
181	<	if ((spt[0] < 0.1) | (spt[0] >= 0.9))
182	<	spt[0] = 0.1 + 0.8*frandom();
183	<	if ((spt[1] < 0.1) | (spt[1] >= 0.9))
184	<	spt[1] = 0.1 + 0.8*frandom();
185	<	}
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);
126	<	dimlist[ndims++] = ambndx(hp,i,j) + 90171;
127	<	rayvalue(arp);                  /* evaluate ray */
128	<	ndims--;                        /* apply coefficient */
129	<	multcolor(arp->rcol, arp->rcoef);
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(&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
201	–	static AMBSAMP *
202	–	ambsample(                              /* initial ambient division sample */
203	–	AMBHEMI *hp,
204	–	int     i,
205	–	int     j
206	–	)
207	–	{
208	–	AMBSAMP *ap = &ambsam(hp,i,j);
209	–	RAY     ar;
210	–	/* generate hemispherical sample */
211	–	if (!getambsamp(&ar, hp, i, j, 0) || ar.rt <= FTINY) {
212	–	memset(ap, 0, sizeof(AMBSAMP));
213	–	return(NULL);
214	–	}
215	–	ap->d = 1.0/ar.rt;              /* limit vertex distance */
216	–	if (ar.rt > 10.0*thescene.cusize)
217	–	ar.rt = 10.0*thescene.cusize;
218	–	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
219	–	copycolor(ap->v, ar.rcol);
220	–	return(ap);
221	–	}
222	–
223	–
159		/* Estimate errors based on 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;
#	Line 238 | Line 174 | getambdiffs(AMBHEMI *hp)
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);
177	>	d2 = normf*(b - bright(ap[-hp->ns].v));
178		d2 *= d2;
179		ep[0] += d2;
180		ep[-hp->ns] += d2;
181		}
182	<	if (j) {                /* from behind */
183	<	d2 = b - bright(ap[-1].v);
184	<	d2 *= d2;
185	<	ep[0] += d2;
186	<	ep[-1] += d2;
187	<	}
182	>	if (!j) continue;
183	>	/* from behind */
184	>	d2 = normf*(b - bright(ap[-1].v));
185	>	d2 *= d2;
186	>	ep[0] += d2;
187	>	ep[-1] += d2;
188	>	if (!i) continue;
189	>	/* diagonal */
190	>	d2 = normf*(b - bright(ap[-hp->ns-1].v));
191	>	d2 *= d2;
192	>	ep[0] += d2;
193	>	ep[-hp->ns-1] += d2;
194		}
195		/* correct for number of neighbors */
196	<	earr[0] *= 2.f;
197	<	earr[hp->ns-1] *= 2.f;
198	<	earr[(hp->ns-1)*hp->ns] *= 2.f;
199	<	earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f;
196	>	earr[0] *= 8./3.;
197	>	earr[hp->ns-1] *= 8./3.;
198	>	earr[(hp->ns-1)*hp->ns] *= 8./3.;
199	>	earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 8./3.;
200		for (i = 1; i < hp->ns-1; i++) {
201	<	earr[i*hp->ns] *= 4./3.;
202	<	earr[i*hp->ns + hp->ns-1] *= 4./3.;
201	>	earr[i*hp->ns] *= 8./5.;
202	>	earr[i*hp->ns + hp->ns-1] *= 8./5.;
203		}
204		for (j = 1; j < hp->ns-1; j++) {
205	<	earr[j] *= 4./3.;
206	<	earr[(hp->ns-1)*hp->ns + j] *= 4./3.;
205	>	earr[j] *= 8./5.;
206	>	earr[(hp->ns-1)*hp->ns + j] *= 8./5.;
207		}
208		return(earr);
209		}
#	Line 269 | Line 211 | getambdiffs(AMBHEMI *hp)
211
212		/* Perform super-sampling on hemisphere (introduces bias) */
213		static void
214	<	ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
214	>	ambsupersamp(AMBHEMI *hp, int cnt)
215		{
216		float   *earr = getambdiffs(hp);
217	<	double  e2sum = 0.0;
217	>	double  e2rem = 0;
218		AMBSAMP *ap;
277	–	RAY     ar;
278	–	double  asum[3];
219		float   *ep;
220	<	int     i, j, n;
220	>	int     i, j, n, nss;
221
222		if (earr == NULL)               /* just skip calc. if no memory */
223		return;
224	<	/* add up estimated variances */
225	<	for (ep = earr + hp->ns*hp->ns; ep-- > earr; )
226	<	e2sum += *ep;
224	>	/* accumulate estimated variances */
225	>	for (ep = earr + hp->ns*hp->ns; ep > earr; )
226	>	e2rem += *--ep;
227		ep = earr;                      /* perform super-sampling */
228		for (ap = hp->sa, i = 0; i < hp->ns; i++)
229		for (j = 0; j < hp->ns; j++, ap++) {
230	<	int     nss = *ep/e2sum*cnt + frandom();
231	<	asum[0] = asum[1] = asum[2] = 0.0;
232	<	for (n = 1; n <= nss; n++) {
233	<	if (!getambsamp(&ar, hp, i, j, n)) {
234	<	nss = n-1;
235	<	break;
296	<	}
297	<	addcolor(asum, ar.rcol);
298	<	}
299	<	if (nss) {              /* update returned ambient value */
300	<	const double    ssf = 1./(nss + 1);
301	<	for (n = 3; n--; )
302	<	acol[n] += ssf*asum[n] +
303	<	(ssf - 1.)*colval(ap->v,n);
304	<	}
305	<	e2sum -= *ep++;         /* update remainders */
306	<	cnt -= nss;
230	>	if (e2rem <= FTINY)
231	>	goto done;      /* nothing left to do */
232	>	nss = *ep/e2rem*cnt + frandom();
233	>	for (n = 1; n <= nss && ambsample(hp,i,j,n); n++)
234	>	if (!--cnt) goto done;
235	>	e2rem -= *ep++;         /* update remainder */
236		}
237	+	done:
238		free(earr);
239		}
240
241
242	<	/* Compute vertex flags, indicating farthest in each direction */
243	<	static uby8 *
244	<	vertex_flags(AMBHEMI *hp)
242	>	static AMBHEMI *
243	>	samp_hemi(                              /* sample indirect hemisphere */
244	>	COLOR   rcol,
245	>	RAY     *r,
246	>	double  wt
247	>	)
248		{
249	<	uby8    *vflags = (uby8 *)calloc(hp->ns*hp->ns, sizeof(uby8));
250	<	uby8    *vf;
251	<	AMBSAMP *ap;
252	<	int     i, j;
253	<
254	<	if (vflags == NULL)
255	<	error(SYSTEM, "out of memory in vertex_flags()");
256	<	vf = vflags;
257	<	ap = hp->sa;            /* compute farthest along first row */
258	<	for (j = 0; j < hp->ns-1; j++, vf++, ap++)
259	<	if (ap[0].d <= ap[1].d)
260	<	vf[0] |= 1<<VDB_X;
261	<	else
262	<	vf[1] |= 1<<VDB_x;
263	<	++vf; ++ap;
264	<	/* flag subsequent rows */
265	<	for (i = 1; i < hp->ns; i++) {
266	<	for (j = 0; j < hp->ns-1; j++, vf++, ap++) {
267	<	if (ap[0].d <= ap[-hp->ns].d)   /* row before */
268	<	vf[0] |= 1<<VDB_y;
269	<	else
270	<	vf[-hp->ns] |= 1<<VDB_Y;
271	<	if (ap[0].d <= ap[1-hp->ns].d)  /* diagonal we care about */
272	<	vf[0] |= 1<<VDB_Xy;
273	<	else
274	<	vf[1-hp->ns] |= 1<<VDB_xY;
275	<	if (ap[0].d <= ap[1].d)         /* column after */
276	<	vf[0] |= 1<<VDB_X;
277	<	else
278	<	vf[1] |= 1<<VDB_x;
279	<	}
280	<	if (ap[0].d <= ap[-hp->ns].d)       /* final column edge */
281	<	vf[0] |= 1<<VDB_y;
282	<	else
283	<	vf[-hp->ns] |= 1<<VDB_Y;
284	<	++vf; ++ap;
249	>	AMBHEMI *hp;
250	>	double  d;
251	>	int     n, i, j;
252	>	/* insignificance check */
253	>	if (bright(rcol) <= FTINY)
254	>	return(NULL);
255	>	/* set number of divisions */
256	>	if (ambacc <= FTINY &&
257	>	wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight)))
258	>	wt = d;                 /* avoid ray termination */
259	>	n = sqrt(ambdiv * wt) + 0.5;
260	>	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
261	>	if (n < i)
262	>	n = i;
263	>	/* allocate sampling array */
264	>	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
265	>	if (hp == NULL)
266	>	error(SYSTEM, "out of memory in samp_hemi");
267	>	hp->rp = r;
268	>	hp->ns = n;
269	>	hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0;
270	>	memset(hp->sa, 0, sizeof(AMBSAMP)*n*n);
271	>	hp->sampOK = 0;
272	>	/* assign coefficient */
273	>	copycolor(hp->acoef, rcol);
274	>	d = 1.0/(n*n);
275	>	scalecolor(hp->acoef, d);
276	>	/* make tangent plane axes */
277	>	if (!getperpendicular(hp->ux, r->ron, 1))
278	>	error(CONSISTENCY, "bad ray direction in samp_hemi");
279	>	VCROSS(hp->uy, r->ron, hp->ux);
280	>	/* sample divisions */
281	>	for (i = hp->ns; i--; )
282	>	for (j = hp->ns; j--; )
283	>	hp->sampOK += ambsample(hp, i, j, 0);
284	>	copycolor(rcol, hp->acol);
285	>	if (!hp->sampOK) {              /* utter failure? */
286	>	free(hp);
287	>	return(NULL);
288		}
289	<	return(vflags);
289	>	if (hp->sampOK < hp->ns*hp->ns) {
290	>	hp->sampOK *= -1;       /* soft failure */
291	>	return(hp);
292	>	}
293	>	n = ambssamp*wt + 0.5;
294	>	if (n > 8) {                    /* perform super-sampling? */
295	>	ambsupersamp(hp, n);
296	>	copycolor(rcol, hp->acol);
297	>	}
298	>	return(hp);                     /* all is well */
299		}
300
301
302		/* Return brightness of farthest ambient sample */
303		static double
304	<	back_ambval(AMBHEMI *hp, int i, int j, int dbit1, int dbit2, const uby8 *vflags)
304	>	back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3)
305		{
306	<	const int       v0 = ambndx(hp,i,j);
307	<	const int       tflags = (1<<dbit1 | 1<<dbit2);
308	<	int             v1, v2;
309	<
310	<	if ((vflags[v0] & tflags) == tflags)    /* is v0 the farthest? */
311	<	return(colval(hp->sa[v0].v,CIEY));
312	<	v1 = adjacent_verti(hp, i, j, dbit1);
313	<	if (vflags[v0] & 1<<dbit2)              /* v1 farthest if v0>v2 */
369	<	return(colval(hp->sa[v1].v,CIEY));
370	<	v2 = adjacent_verti(hp, i, j, dbit2);
371	<	if (vflags[v0] & 1<<dbit1)              /* v2 farthest if v0>v1 */
372	<	return(colval(hp->sa[v2].v,CIEY));
373	<	/* else check if v1>v2 */
374	<	if (vflags[v1] & 1<<vdb_edge(dbit1,dbit2))
375	<	return(colval(hp->sa[v1].v,CIEY));
376	<	return(colval(hp->sa[v2].v,CIEY));
306	>	if (hp->sa[n1].d <= hp->sa[n2].d) {
307	>	if (hp->sa[n1].d <= hp->sa[n3].d)
308	>	return(colval(hp->sa[n1].v,CIEY));
309	>	return(colval(hp->sa[n3].v,CIEY));
310	>	}
311	>	if (hp->sa[n2].d <= hp->sa[n3].d)
312	>	return(colval(hp->sa[n2].v,CIEY));
313	>	return(colval(hp->sa[n3].v,CIEY));
314		}
315
316
317		/* Compute vectors and coefficients for Hessian/gradient calcs */
318		static void
319	<	comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags)
319	>	comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1)
320		{
321	<	const int       i0 = ambndx(hp,i,j);
322	<	double          rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
386	<	int             i1, ii;
321	>	double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
322	>	int     ii;
323
324	<	ftp->valid = 0;                 /* check if we can skip this edge */
325	<	ii = adjacent_trifl[dbit];
326	<	if ((vflags[i0] & ii) == ii)    /* cancels if vertex used as value */
391	<	return;
392	<	i1 = adjacent_verti(hp, i, j, dbit);
393	<	ii = adjacent_trifl[VDB_OPP(dbit)];
394	<	if ((vflags[i1] & ii) == ii)    /* on either end (for both triangles) */
395	<	return;
396	<	/* else go ahead with calculation */
397	<	VSUB(ftp->r_i, hp->sa[i0].p, hp->rp->rop);
398	<	VSUB(ftp->r_i1, hp->sa[i1].p, hp->rp->rop);
399	<	VSUB(ftp->e_i, hp->sa[i1].p, hp->sa[i0].p);
324	>	VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop);
325	>	VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop);
326	>	VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p);
327		VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
328		rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
329		dot_e = DOT(ftp->e_i,ftp->e_i);
#	Line 410 | Line 337 | comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int
337		J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
338		for (ii = 3; ii--; )
339		ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
413	–	ftp->valid++;
340		}
341
342
#	Line 436 | Line 362 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
362		double  d1, d2, d3, d4;
363		double  I3, J3, K3;
364		int     i, j;
439	–
440	–	if (!ftp->valid) {              /* preemptive test */
441	–	memset(hess, 0, sizeof(FVECT)*3);
442	–	return;
443	–	}
365		/* compute intermediate coefficients */
366		d1 = 1.0/DOT(ftp->r_i,ftp->r_i);
367		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
#	Line 504 | Line 425 | comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
425		double  f1;
426		int     i;
427
507	–	if (!ftp->valid) {              /* preemptive test */
508	–	memset(grad, 0, sizeof(FVECT));
509	–	return;
510	–	}
428		f1 = 2.0*DOT(nrm, ftp->rcp);
429		VCROSS(ncp, nrm, ftp->e_i);
430		for (i = 3; i--; )
#	Line 537 | Line 454 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
454
455
456		/* Compute anisotropic radii and eigenvector directions */
457	<	static int
457	>	static void
458		eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
459		{
460		double  hess2[2][2];
#	Line 559 | Line 476 | eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
476		if (i == 1)                     /* double-root (circle) */
477		evalue[1] = evalue[0];
478		if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) |
479	<	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) )
480	<	error(INTERNAL, "bad eigenvalue calculation");
481	<
479	>	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) {
480	>	ra[0] = ra[1] = maxarad;
481	>	return;
482	>	}
483		if (evalue[0] > evalue[1]) {
484		ra[0] = sqrt(sqrt(4.0/evalue[0]));
485		ra[1] = sqrt(sqrt(4.0/evalue[1]));
#	Line 596 | Line 514 | ambHessian(                            /* anisotropic radii & pos. gradient */
514		static char     memerrmsg[] = "out of memory in ambHessian()";
515		FVECT           (*hessrow)[3] = NULL;
516		FVECT           *gradrow = NULL;
599	–	uby8            *vflags;
517		FVECT           hessian[3];
518		FVECT           gradient;
519		FFTRI           fftr;
#	Line 618 | Line 535 | ambHessian(                            /* anisotropic radii & pos. gradient */
535		error(SYSTEM, memerrmsg);
536		memset(gradient, 0, sizeof(gradient));
537		}
621	–	/* get vertex position flags */
622	–	vflags = vertex_flags(hp);
538		/* compute first row of edges */
539		for (j = 0; j < hp->ns-1; j++) {
540	<	comp_fftri(&fftr, hp, 0, j, VDB_X, vflags);
540	>	comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1));
541		if (hessrow != NULL)
542		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
543		if (gradrow != NULL)
#	Line 632 | Line 547 | ambHessian(                            /* anisotropic radii & pos. gradient */
547		for (i = 0; i < hp->ns-1; i++) {
548		FVECT       hesscol[3];     /* compute first vertical edge */
549		FVECT       gradcol;
550	<	comp_fftri(&fftr, hp, i, 0, VDB_Y, vflags);
550	>	comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0));
551		if (hessrow != NULL)
552		comp_hessian(hesscol, &fftr, hp->rp->ron);
553		if (gradrow != NULL)
#	Line 641 | Line 556 | ambHessian(                            /* anisotropic radii & pos. gradient */
556		FVECT   hessdia[3];     /* compute triangle contributions */
557		FVECT   graddia;
558		double  backg;
559	<	backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags);
559	>	backg = back_ambval(hp, AI(hp,i,j),
560	>	AI(hp,i,j+1), AI(hp,i+1,j));
561		/* diagonal (inner) edge */
562	<	comp_fftri(&fftr, hp, i, j+1, VDB_xY, vflags);
562	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j));
563		if (hessrow != NULL) {
564		comp_hessian(hessdia, &fftr, hp->rp->ron);
565		rev_hessian(hesscol);
#	Line 655 | Line 571 | ambHessian(                            /* anisotropic radii & pos. gradient */
571		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
572		}
573		/* initialize edge in next row */
574	<	comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags);
574	>	comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j));
575		if (hessrow != NULL)
576		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
577		if (gradrow != NULL)
578		comp_gradient(gradrow[j], &fftr, hp->rp->ron);
579		/* new column edge & paired triangle */
580	<	backg = back_ambval(hp, i+1, j+1, VDB_x, VDB_y, vflags);
581	<	comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags);
580	>	backg = back_ambval(hp, AI(hp,i+1,j+1),
581	>	AI(hp,i+1,j), AI(hp,i,j+1));
582	>	comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1));
583		if (hessrow != NULL) {
584		comp_hessian(hesscol, &fftr, hp->rp->ron);
585		rev_hessian(hessdia);
#	Line 682 | Line 599 | ambHessian(                            /* anisotropic radii & pos. gradient */
599		/* release row buffers */
600		if (hessrow != NULL) free(hessrow);
601		if (gradrow != NULL) free(gradrow);
685	–	free(vflags);
602
603		if (ra != NULL)                 /* extract eigenvectors & radii */
604		eigenvectors(uv, ra, hessian);
#	Line 723 | Line 639 | static uint32
639		ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1)
640		{
641		const double    max_d = 1.0/(minarad*ambacc + 0.001);
642	<	const double    ang_res = 0.5*PI/(hp->ns-1);
643	<	const double    ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY));
642	>	const double    ang_res = 0.5*PI/hp->ns;
643	>	const double    ang_step = ang_res/((int)(16/PI*ang_res) + 1.01);
644		double          avg_d = 0;
645		uint32          flgs = 0;
646	+	FVECT           vec;
647	+	double          u, v;
648	+	double          ang, a1;
649	+	OBJREC          *m;
650		int             i, j;
651	<	/* check distances above us */
651	>	/* don't bother for a few samples */
652	>	if (hp->ns < 8)
653	>	return(0);
654	>	/* check distances overhead */
655		for (i = hp->ns*3/4; i-- > hp->ns>>2; )
656		for (j = hp->ns*3/4; j-- > hp->ns>>2; )
657		avg_d += ambsam(hp,i,j).d;
658		avg_d *= 4.0/(hp->ns*hp->ns);
659	<	if (avg_d >= max_d)             /* too close to corral? */
659	>	if (avg_d*r0 >= 1.0)            /* ceiling too low for corral? */
660		return(0);
661	+	if (avg_d >= max_d)             /* insurance */
662	+	return(0);
663		/* else circle around perimeter */
664		for (i = 0; i < hp->ns; i++)
665		for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) {
666		AMBSAMP *ap = &ambsam(hp,i,j);
742	–	FVECT   vec;
743	–	double  u, v;
744	–	double  ang, a1;
745	–	int     abp;
667		if ((ap->d <= FTINY) | (ap->d >= max_d))
668		continue;       /* too far or too near */
669		VSUB(vec, ap->p, hp->rp->rop);
670	<	u = DOT(vec, uv[0]) * ap->d;
671	<	v = DOT(vec, uv[1]) * ap->d;
672	<	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0)
670	>	u = DOT(vec, uv[0]);
671	>	v = DOT(vec, uv[1]);
672	>	if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v)
673		continue;       /* occluder outside ellipse */
674		ang = atan2a(v, u);     /* else set direction flags */
675	<	for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step)
675	>	for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step)
676		flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0)));
677		}
678		return(flgs);
#	Line 770 | Line 691 | doambient(                             /* compute ambient component */
691		uint32  *crlp                   /* returned (optional) */
692		)
693		{
694	<	AMBHEMI *hp = inithemi(rcol, r, wt);
774	<	int     cnt;
694	>	AMBHEMI *hp = samp_hemi(rcol, r, wt);
695		FVECT   my_uv[2];
696	<	double  d, K, acol[3];
696	>	double  d, K;
697		AMBSAMP *ap;
698	<	int     i, j;
699	<	/* check/initialize */
780	<	if (hp == NULL)
781	<	return(0);
698	>	int     i;
699	>	/* clear return values */
700		if (uv != NULL)
701		memset(uv, 0, sizeof(FVECT)*2);
702		if (ra != NULL)
#	Line 789 | Line 707 | doambient(                             /* compute ambient component */
707		dg[0] = dg[1] = 0.0;
708		if (crlp != NULL)
709		*crlp = 0;
710	<	/* sample the hemisphere */
711	<	acol[0] = acol[1] = acol[2] = 0.0;
712	<	cnt = 0;
713	<	for (i = hp->ns; i--; )
714	<	for (j = hp->ns; j--; )
715	<	if ((ap = ambsample(hp, i, j)) != NULL) {
716	<	addcolor(acol, ap->v);
799	<	++cnt;
800	<	}
801	<	if (!cnt) {
802	<	setcolor(rcol, 0.0, 0.0, 0.0);
803	<	free(hp);
804	<	return(0);              /* no valid samples */
710	>	if (hp == NULL)                 /* sampling falure? */
711	>	return(0);
712	>
713	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL) ||
714	>	(hp->sampOK < 0) | (hp->ns < 6)) {
715	>	free(hp);               /* Hessian not requested/possible */
716	>	return(-1);             /* value-only return value */
717		}
718	<	if (cnt < hp->ns*hp->ns) {      /* incomplete sampling? */
807	<	copycolor(rcol, acol);
808	<	free(hp);
809	<	return(-1);             /* return value w/o Hessian */
810	<	}
811	<	cnt = ambssamp*wt + 0.5;        /* perform super-sampling? */
812	<	if (cnt > 8)
813	<	ambsupersamp(acol, hp, cnt);
814	<	copycolor(rcol, acol);          /* final indirect irradiance/PI */
815	<	if ((ra == NULL) & (pg == NULL) & (dg == NULL)) {
816	<	free(hp);
817	<	return(-1);             /* no radius or gradient calc. */
818	<	}
819	<	if ((d = bright(acol)) > FTINY) {       /* normalize Y values */
718	>	if ((d = bright(rcol)) > FTINY) {       /* normalize Y values */
719		d = 0.99*(hp->ns*hp->ns)/d;
720		K = 0.01;
721		} else {                        /* or fall back on geometric Hessian */
722		K = 1.0;
723		pg = NULL;
724		dg = NULL;
725	+	crlp = NULL;
726		}
727		ap = hp->sa;                    /* relative Y channel from here on... */
728		for (i = hp->ns*hp->ns; i--; ap++)
#	Line 850 | Line 750 | doambient(                             /* compute ambient component */
750		if (ra[1] < minarad)
751		ra[1] = minarad;
752		}
753	<	ra[0] *= d = 1.0/sqrt(sqrt(wt));
753	>	ra[0] *= d = 1.0/sqrt(wt);
754		if ((ra[1] *= d) > 2.0*ra[0])
755		ra[1] = 2.0*ra[0];
756		if (ra[1] > maxarad) {
#	Line 858 | Line 758 | doambient(                             /* compute ambient component */
758		if (ra[0] > maxarad)
759		ra[0] = maxarad;
760		}
761	<	if (crlp != NULL)       /* flag encroached directions */
761	>	/* flag encroached directions */
762	>	if (crlp != NULL)
763		*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc);
764		if (pg != NULL) {       /* cap gradient if necessary */
765		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)