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.32 by greg, Thu Apr 24 17:36:43 2014 UTC vs.
Revision 2.44 by greg, Thu May 1 16:06:11 2014 UTC

#	Line 22 | Line 22 | static const char      RCSid[] = "$Id$";
22		extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
23
24		typedef struct {
25	+	COLOR   v;              /* hemisphere sample value */
26	+	FVECT   p;              /* intersection point */
27	+	} AMBSAMP;              /* sample value */
28	+
29	+	typedef struct {
30		RAY     *rp;            /* originating ray sample */
31		FVECT   ux, uy;         /* tangent axis unit vectors */
32		int     ns;             /* number of samples per axis */
33		COLOR   acoef;          /* division contribution coefficient */
34	<	struct s_ambsamp {
30	<	COLOR   v;              /* hemisphere sample value */
31	<	FVECT   p;              /* intersection point */
32	<	} sa[1];                /* sample array (extends struct) */
34	>	AMBSAMP sa[1];          /* sample array (extends struct) */
35		}  AMBHEMI;             /* ambient sample hemisphere */
36
37	<	#define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
37	>	#define ambsam(h,i,j)   (h)->sa[(i)*(h)->ns + (j)]
38
39		typedef struct {
40	<	FVECT   r_i, r_i1, e_i, rI2_eJ2;
41	<	double  nf, I1, I2;
40	>	FVECT   r_i, r_i1, e_i, rcp, rI2_eJ2;
41	>	double  I1, I2;
42		} FFTRI;                /* vectors and coefficients for Hessian calculation */
43
44
#	Line 59 | Line 61 | inithemi(                      /* initialize sampling hemisphere */
61		if (n < i)
62		n = i;
63		/* allocate sampling array */
64	<	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
63	<	sizeof(struct s_ambsamp)*(n*n - 1));
64	>	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
65		if (hp == NULL)
66		return(NULL);
67		hp->rp = r;
#	Line 70 | Line 71 | inithemi(                      /* initialize sampling hemisphere */
71		d = 1.0/(n*n);
72		scalecolor(hp->acoef, d);
73		/* make tangent plane axes */
74	<	hp->uy[0] = 0.1 - 0.2*frandom();
75	<	hp->uy[1] = 0.1 - 0.2*frandom();
76	<	hp->uy[2] = 0.1 - 0.2*frandom();
77	<	for (i = 0; i < 3; i++)
78	<	if (r->ron[i] < 0.6 && r->ron[i] > -0.6)
74	>	hp->uy[0] = 0.5 - frandom();
75	>	hp->uy[1] = 0.5 - frandom();
76	>	hp->uy[2] = 0.5 - frandom();
77	>	for (i = 3; i--; )
78	>	if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
79		break;
80	<	if (i >= 3)
81	<	error(CONSISTENCY, "bad ray direction in inithemi()");
80	>	if (i < 0)
81	>	error(CONSISTENCY, "bad ray direction in inithemi");
82		hp->uy[i] = 1.0;
83		VCROSS(hp->ux, hp->uy, r->ron);
84		normalize(hp->ux);
#	Line 87 | Line 88 | inithemi(                      /* initialize sampling hemisphere */
88		}
89
90
91	<	static struct s_ambsamp *
92	<	ambsample(                              /* sample an ambient direction */
93	<	AMBHEMI *hp,
93	<	int     i,
94	<	int     j
95	<	)
91	>	/* Sample ambient division and apply weighting coefficient */
92	>	static int
93	>	getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
94		{
95	<	struct s_ambsamp        *ap = &ambsamp(hp,i,j);
96	<	RAY                     ar;
99	<	double                  spt[2], zd;
100	<	int                     ii;
95	>	int     hlist[3], ii;
96	>	double  spt[2], zd;
97		/* ambient coefficient for weight */
98		if (ambacc > FTINY)
99	<	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
99	>	setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL);
100		else
101	<	copycolor(ar.rcoef, hp->acoef);
102	<	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
103	<	goto badsample;
101	>	copycolor(arp->rcoef, hp->acoef);
102	>	if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0)
103	>	return(0);
104		if (ambacc > FTINY) {
105	<	multcolor(ar.rcoef, hp->acoef);
106	<	scalecolor(ar.rcoef, 1./AVGREFL);
105	>	multcolor(arp->rcoef, hp->acoef);
106	>	scalecolor(arp->rcoef, 1./AVGREFL);
107		}
108	<	/* generate hemispherical sample */
109	<	SDsquare2disk(spt,      (i+.1+.8*frandom())/hp->ns,
110	<	(j+.1+.8*frandom())/hp->ns );
108	>	hlist[0] = hp->rp->rno;
109	>	hlist[1] = i;
110	>	hlist[2] = j;
111	>	multisamp(spt, 2, urand(ilhash(hlist,3)+n));
112	>	if (!n) {                       /* avoid border samples for n==0 */
113	>	if ((spt[0] < 0.1) | (spt[0] > 0.9))
114	>	spt[0] = 0.1 + 0.8*frandom();
115	>	if ((spt[1] < 0.1) | (spt[1] > 0.9))
116	>	spt[1] = 0.1 + 0.8*frandom();
117	>	}
118	>	SDsquare2disk(spt, (i+spt[0])/hp->ns, (j+spt[1])/hp->ns);
119		zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
120		for (ii = 3; ii--; )
121	<	ar.rdir[ii] =   spt[0]*hp->ux[ii] +
121	>	arp->rdir[ii] = spt[0]*hp->ux[ii] +
122		spt[1]*hp->uy[ii] +
123		zd*hp->rp->ron[ii];
124	<	checknorm(ar.rdir);
124	>	checknorm(arp->rdir);
125		dimlist[ndims++] = i*hp->ns + j + 90171;
126	<	rayvalue(&ar);                  /* evaluate ray */
127	<	ndims--;
128	<	if (ar.rt > 20.0*maxarad)       /* limit vertex distance */
129	<	ar.rt = 20.0*maxarad;
126	>	rayvalue(arp);                  /* evaluate ray */
127	>	ndims--;                        /* apply coefficient */
128	>	multcolor(arp->rcol, arp->rcoef);
129	>	return(1);
130	>	}
131	>
132	>
133	>	static AMBSAMP *
134	>	ambsample(                              /* initial ambient division sample */
135	>	AMBHEMI *hp,
136	>	int     i,
137	>	int     j
138	>	)
139	>	{
140	>	AMBSAMP *ap = &ambsam(hp,i,j);
141	>	RAY     ar;
142	>	/* generate hemispherical sample */
143	>	if (!getambsamp(&ar, hp, i, j, 0))
144	>	goto badsample;
145	>	/* limit vertex distance */
146	>	if (ar.rt > 10.0*thescene.cusize)
147	>	ar.rt = 10.0*thescene.cusize;
148		else if (ar.rt <= FTINY)        /* should never happen! */
149		goto badsample;
150		VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
129	–	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
151		copycolor(ap->v, ar.rcol);
152		return(ap);
153		badsample:
#	Line 136 | Line 157 | badsample:
157		}
158
159
160	+	/* Estimate errors based on ambient division differences */
161	+	static float *
162	+	getambdiffs(AMBHEMI *hp)
163	+	{
164	+	float   *earr = calloc(hp->ns*hp->ns, sizeof(float));
165	+	float   *ep;
166	+	AMBSAMP *ap;
167	+	double  b, d2;
168	+	int     i, j;
169	+
170	+	if (earr == NULL)               /* out of memory? */
171	+	return(NULL);
172	+	/* compute 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;
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	+	}
188	+	}
189	+	/* correct for number of neighbors */
190	+	earr[0] *= 2.f;
191	+	earr[hp->ns-1] *= 2.f;
192	+	earr[(hp->ns-1)*hp->ns] *= 2.f;
193	+	earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f;
194	+	for (i = 1; i < hp->ns-1; i++) {
195	+	earr[i*hp->ns] *= 4./3.;
196	+	earr[i*hp->ns + hp->ns-1] *= 4./3.;
197	+	}
198	+	for (j = 1; j < hp->ns-1; j++) {
199	+	earr[j] *= 4./3.;
200	+	earr[(hp->ns-1)*hp->ns + j] *= 4./3.;
201	+	}
202	+	return(earr);
203	+	}
204	+
205	+
206	+	/* Perform super-sampling on hemisphere (introduces bias) */
207	+	static void
208	+	ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
209	+	{
210	+	float   *earr = getambdiffs(hp);
211	+	double  e2sum = 0;
212	+	AMBSAMP *ap;
213	+	RAY     ar;
214	+	COLOR   asum;
215	+	float   *ep;
216	+	int     i, j, n;
217	+
218	+	if (earr == NULL)               /* just skip calc. if no memory */
219	+	return;
220	+	/* add up estimated variances */
221	+	for (ep = earr + hp->ns*hp->ns; ep-- > earr; )
222	+	e2sum += *ep;
223	+	ep = earr;                      /* perform super-sampling */
224	+	for (ap = hp->sa, i = 0; i < hp->ns; i++)
225	+	for (j = 0; j < hp->ns; j++, ap++) {
226	+	int     nss = *ep/e2sum*cnt + frandom();
227	+	setcolor(asum, 0., 0., 0.);
228	+	for (n = 1; n <= nss; n++) {
229	+	if (!getambsamp(&ar, hp, i, j, n)) {
230	+	nss = n-1;
231	+	break;
232	+	}
233	+	addcolor(asum, ar.rcol);
234	+	}
235	+	if (nss) {              /* update returned ambient value */
236	+	const double    ssf = 1./(nss + 1);
237	+	for (n = 3; n--; )
238	+	acol[n] += ssf*colval(asum,n) +
239	+	(ssf - 1.)*colval(ap->v,n);
240	+	}
241	+	e2sum -= *ep++;         /* update remainders */
242	+	cnt -= nss;
243	+	}
244	+	free(earr);
245	+	}
246	+
247	+
248		/* Compute vectors and coefficients for Hessian/gradient calcs */
249		static void
250		comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop)
251		{
252	<	FVECT   vcp;
144	<	double  dot_e, dot_er, rdot_r, rdot_r1, J2;
252	>	double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
253		int     i;
254
255		VSUB(ftp->r_i, ap0, rop);
256		VSUB(ftp->r_i1, ap1, rop);
257		VSUB(ftp->e_i, ap1, ap0);
258	<	VCROSS(vcp, ftp->e_i, ftp->r_i);
259	<	ftp->nf = 1.0/DOT(vcp,vcp);
258	>	VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
259	>	rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
260		dot_e = DOT(ftp->e_i,ftp->e_i);
261		dot_er = DOT(ftp->e_i, ftp->r_i);
262		rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i);
263		rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
264		ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) *
265	<	sqrt( ftp->nf );
265	>	sqrt( rdot_cp );
266		ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
267	<	dot_e*ftp->I1 )*0.5*ftp->nf;
267	>	dot_e*ftp->I1 )*0.5*rdot_cp;
268		J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
269		for (i = 3; i--; )
270		ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i];
#	Line 180 | Line 288 | compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
288		static void
289		comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
290		{
291	<	FVECT   vcp;
291	>	FVECT   ncp;
292		FVECT   m1[3], m2[3], m3[3], m4[3];
293		double  d1, d2, d3, d4;
294		double  I3, J3, K3;
#	Line 190 | Line 298 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
298		d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
299		d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
300		d4 = DOT(ftp->e_i, ftp->r_i);
301	<	I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 +
302	<	3.0/d3*ftp->I2 );
301	>	I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 )
302	>	/ ( 4.0*DOT(ftp->rcp,ftp->rcp) );
303		J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3;
304		K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3);
305		/* intermediate matrices */
306	<	VCROSS(vcp, nrm, ftp->e_i);
307	<	compose_matrix(m1, vcp, ftp->rI2_eJ2);
306	>	VCROSS(ncp, nrm, ftp->e_i);
307	>	compose_matrix(m1, ncp, ftp->rI2_eJ2);
308		compose_matrix(m2, ftp->r_i, ftp->r_i);
309		compose_matrix(m3, ftp->e_i, ftp->e_i);
310		compose_matrix(m4, ftp->r_i, ftp->e_i);
311	<	VCROSS(vcp, ftp->r_i, ftp->e_i);
204	<	d1 = DOT(nrm, vcp);
311	>	d1 = DOT(nrm, ftp->rcp);
312		d2 = -d1*ftp->I2;
313		d1 *= 2.0;
314		for (i = 3; i--; )              /* final matrix sum */
#	Line 245 | Line 352 | add2hessian(FVECT hess[3], FVECT ehess1[3],
352		static void
353		comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
354		{
355	<	FVECT   vcp;
355	>	FVECT   ncp;
356		double  f1;
357		int     i;
358
359	<	VCROSS(vcp, ftp->r_i, ftp->r_i1);
360	<	f1 = 2.0*DOT(nrm, vcp);
254	<	VCROSS(vcp, nrm, ftp->e_i);
359	>	f1 = 2.0*DOT(nrm, ftp->rcp);
360	>	VCROSS(ncp, nrm, ftp->e_i);
361		for (i = 3; i--; )
362	<	grad[i] = (-0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] );
362	>	grad[i] = (-0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
363		}
364
365
#	Line 280 | Line 386 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
386
387		/* Return brightness of furthest ambient sample */
388		static COLORV
389	<	back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2,
284	<	struct s_ambsamp *ap3, FVECT orig)
389	>	back_ambval(AMBSAMP *ap1, AMBSAMP *ap2, AMBSAMP *ap3, FVECT orig)
390		{
391		COLORV  vback;
392		FVECT   vec;
#	Line 321 | Line 426 | eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
426		hess2[0][1] = DOT(uv[0], b);
427		hess2[1][0] = DOT(uv[1], a);
428		hess2[1][1] = DOT(uv[1], b);
429	<	/* compute eigenvalues */
430	<	if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
431	<	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 ||
432	<	(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY ||
433	<	(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY )
429	>	/* compute eigenvalue(s) */
430	>	i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
431	>	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]);
432	>	if (i == 1)                     /* double-root (circle) */
433	>	evalue[1] = evalue[0];
434	>	if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) |
435	>	((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) )
436		error(INTERNAL, "bad eigenvalue calculation");
437
438		if (evalue[0] > evalue[1]) {
#	Line 385 | Line 492 | ambHessian(                            /* anisotropic radii & pos. gradient */
492		}
493		/* compute first row of edges */
494		for (j = 0; j < hp->ns-1; j++) {
495	<	comp_fftri(&fftr, ambsamp(hp,0,j).p,
496	<	ambsamp(hp,0,j+1).p, hp->rp->rop);
495	>	comp_fftri(&fftr, ambsam(hp,0,j).p,
496	>	ambsam(hp,0,j+1).p, hp->rp->rop);
497		if (hessrow != NULL)
498		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
499		if (gradrow != NULL)
#	Line 396 | Line 503 | ambHessian(                            /* anisotropic radii & pos. gradient */
503		for (i = 0; i < hp->ns-1; i++) {
504		FVECT       hesscol[3];     /* compute first vertical edge */
505		FVECT       gradcol;
506	<	comp_fftri(&fftr, ambsamp(hp,i,0).p,
507	<	ambsamp(hp,i+1,0).p, hp->rp->rop);
506	>	comp_fftri(&fftr, ambsam(hp,i,0).p,
507	>	ambsam(hp,i+1,0).p, hp->rp->rop);
508		if (hessrow != NULL)
509		comp_hessian(hesscol, &fftr, hp->rp->ron);
510		if (gradrow != NULL)
#	Line 406 | Line 513 | ambHessian(                            /* anisotropic radii & pos. gradient */
513		FVECT   hessdia[3];     /* compute triangle contributions */
514		FVECT   graddia;
515		COLORV  backg;
516	<	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
517	<	&ambsamp(hp,i+1,j), hp->rp->rop);
516	>	backg = back_ambval(&ambsam(hp,i,j), &ambsam(hp,i,j+1),
517	>	&ambsam(hp,i+1,j), hp->rp->rop);
518		/* diagonal (inner) edge */
519	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
520	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
519	>	comp_fftri(&fftr, ambsam(hp,i,j+1).p,
520	>	ambsam(hp,i+1,j).p, hp->rp->rop);
521		if (hessrow != NULL) {
522		comp_hessian(hessdia, &fftr, hp->rp->ron);
523		rev_hessian(hesscol);
524		add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
525		}
526	<	if (gradient != NULL) {
526	>	if (gradrow != NULL) {
527		comp_gradient(graddia, &fftr, hp->rp->ron);
528		rev_gradient(gradcol);
529		add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
530		}
531		/* initialize edge in next row */
532	<	comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
533	<	ambsamp(hp,i+1,j).p, hp->rp->rop);
532	>	comp_fftri(&fftr, ambsam(hp,i+1,j+1).p,
533	>	ambsam(hp,i+1,j).p, hp->rp->rop);
534		if (hessrow != NULL)
535		comp_hessian(hessrow[j], &fftr, hp->rp->ron);
536		if (gradrow != NULL)
537		comp_gradient(gradrow[j], &fftr, hp->rp->ron);
538		/* new column edge & paired triangle */
539	<	backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
540	<	&ambsamp(hp,i+1,j), hp->rp->rop);
541	<	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
539	>	backg = back_ambval(&ambsam(hp,i,j+1), &ambsam(hp,i+1,j+1),
540	>	&ambsam(hp,i+1,j), hp->rp->rop);
541	>	comp_fftri(&fftr, ambsam(hp,i,j+1).p, ambsam(hp,i+1,j+1).p,
542		hp->rp->rop);
543		if (hessrow != NULL) {
544		comp_hessian(hesscol, &fftr, hp->rp->ron);
#	Line 466 | Line 573 | ambHessian(                            /* anisotropic radii & pos. gradient */
573		static void
574		ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
575		{
576	<	struct s_ambsamp        *ap;
577	<	double                  dgsum[2];
578	<	int                     n;
579	<	FVECT                   vd;
580	<	double                  gfact;
576	>	AMBSAMP *ap;
577	>	double  dgsum[2];
578	>	int     n;
579	>	FVECT   vd;
580	>	double  gfact;
581
582		dgsum[0] = dgsum[1] = 0.0;      /* sum values times -tan(theta) */
583		for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
#	Line 478 | Line 585 | ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
585		VSUB(vd, ap->p, hp->rp->rop);
586		/* brightness over cosine factor */
587		gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
588	<	/* -sine = -proj_radius/vd_length */
589	<	dgsum[0] += DOT(uv[1], vd) * gfact;
590	<	dgsum[1] -= DOT(uv[0], vd) * gfact;
588	>	/* sine = proj_radius/vd_length */
589	>	dgsum[0] -= DOT(uv[1], vd) * gfact;
590	>	dgsum[1] += DOT(uv[0], vd) * gfact;
591		}
592		dg[0] = dgsum[0] / (hp->ns*hp->ns);
593		dg[1] = dgsum[1] / (hp->ns*hp->ns);
#	Line 498 | Line 605 | doambient(                             /* compute ambient component */
605		float   dg[2]                   /* returned (optional) */
606		)
607		{
608	<	AMBHEMI                 *hp = inithemi(rcol, r, wt);
609	<	int                     cnt = 0;
610	<	FVECT                   my_uv[2];
611	<	double                  d, acol[3];
612	<	struct s_ambsamp        *ap;
613	<	int                     i, j;
608	>	AMBHEMI *hp = inithemi(rcol, r, wt);
609	>	int     cnt = 0;
610	>	FVECT   my_uv[2];
611	>	double  d, K, acol[3];
612	>	AMBSAMP *ap;
613	>	int     i, j;
614		/* check/initialize */
615		if (hp == NULL)
616		return(0);
#	Line 528 | Line 635 | doambient(                             /* compute ambient component */
635		free(hp);
636		return(0);              /* no valid samples */
637		}
638	+	if (cnt < hp->ns*hp->ns) {      /* incomplete sampling? */
639	+	copycolor(rcol, acol);
640	+	free(hp);
641	+	return(-1);             /* return value w/o Hessian */
642	+	}
643	+	cnt = ambssamp*wt + 0.5;        /* perform super-sampling? */
644	+	if (cnt > 0)
645	+	ambsupersamp(acol, hp, cnt);
646		copycolor(rcol, acol);          /* final indirect irradiance/PI */
647	<	if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
533	<	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
647	>	if ((ra == NULL) & (pg == NULL) & (dg == NULL)) {
648		free(hp);
649		return(-1);             /* no radius or gradient calc. */
650		}
651	<	if (bright(acol) > FTINY)       /* normalize Y values */
652	<	d = cnt/bright(acol);
653	<	else
651	>	if (bright(acol) > FTINY) {     /* normalize Y values */
652	>	d = 0.99*cnt/bright(acol);
653	>	K = 0.01;
654	>	} else {                        /* geometric Hessian fall-back */
655		d = 0.0;
656	+	K = 1.0;
657	+	pg = NULL;
658	+	dg = NULL;
659	+	}
660		ap = hp->sa;                    /* relative Y channel from here on... */
661		for (i = hp->ns*hp->ns; i--; ap++)
662	<	colval(ap->v,CIEY) = bright(ap->v)*d + 0.01;
662	>	colval(ap->v,CIEY) = bright(ap->v)*d + K;
663
664		if (uv == NULL)                 /* make sure we have axis pointers */
665		uv = my_uv;
#	Line 551 | Line 670 | doambient(                             /* compute ambient component */
670		ambdirgrad(hp, uv, dg);
671
672		if (ra != NULL) {               /* scale/clamp radii */
673	+	if (pg != NULL) {
674	+	if (ra[0]*(d = fabs(pg[0])) > 1.0)
675	+	ra[0] = 1.0/d;
676	+	if (ra[1]*(d = fabs(pg[1])) > 1.0)
677	+	ra[1] = 1.0/d;
678	+	if (ra[0] > ra[1])
679	+	ra[0] = ra[1];
680	+	}
681		if (ra[0] < minarad) {
682		ra[0] = minarad;
683		if (ra[1] < minarad)
#	Line 563 | Line 690 | doambient(                             /* compute ambient component */
690		ra[1] = maxarad;
691		if (ra[0] > maxarad)
692		ra[0] = maxarad;
693	+	}
694	+	if (pg != NULL) {       /* cap gradient if necessary */
695	+	d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];
696	+	if (d > 1.0) {
697	+	d = 1.0/sqrt(d);
698	+	pg[0] *= d;
699	+	pg[1] *= d;
700	+	}
701		}
702		}
703		free(hp);                       /* clean up and return */

Diff Legend

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