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.35 by greg, Fri Apr 25 18:39:22 2014 UTC vs.
Revision 2.43 by greg, Thu May 1 16:01:59 2014 UTC

#	Line 32 | Line 32 | typedef struct {
32		} sa[1];                /* sample array (extends struct) */
33		}  AMBHEMI;             /* ambient sample hemisphere */
34
35	<	#define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
35	>	typedef struct s_ambsamp        AMBSAMP;
36
37	+	#define ambsam(h,i,j)   (h)->sa[(i)*(h)->ns + (j)]
38	+
39		typedef struct {
40		FVECT   r_i, r_i1, e_i, rcp, rI2_eJ2;
41		double  I1, I2;
#	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--;
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);
130	–	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
151		copycolor(ap->v, ar.rcol);
152		return(ap);
153		badsample:
#	Line 137 | 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)
#	Line 278 | 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,
282	<	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 319 | 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) |
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
#	Line 383 | 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 394 | 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 404 | 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 464 | 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 476 | 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 496 | 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 526 | 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? */
531	<	(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;

Diff Legend

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