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 1.4 by greg, Tue Jun 11 16:37:49 1991 UTC vs.
Revision 2.28 by greg, Sat Apr 19 19:20:47 2014 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1991 Regents of the University of California */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ LBL";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		* Routines to compute "ambient" values using Monte Carlo
6	+	*
7	+	*  Hessian calculations based on "Practical Hessian-Based Error Control
8	+	*      for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz
9	+	*      from ACM SIGGRAPH Asia 2012 conference proceedings.
10	+	*
11	+	*  Declarations of external symbols in ambient.h
12		*/
13
14	<	#include  "ray.h"
14	>	#include "copyright.h"
15
16	+	#include  "ray.h"
17		#include  "ambient.h"
14	–
18		#include  "random.h"
19
20	<	typedef struct {
18	<	short  t, p;            /* theta, phi indices */
19	<	COLOR  v;               /* value sum */
20	<	float  r;               /* 1/distance sum */
21	<	float  k;               /* variance for this division */
22	<	int  n;                 /* number of subsamples */
23	<	}  AMBSAMP;             /* ambient sample division */
20	>	#ifdef NEWAMB
21
22	+	extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
23	+
24		typedef struct {
25	<	FVECT  ux, uy, uz;      /* x, y and z axis directions */
26	<	short  nt, np;          /* number of theta and phi directions */
25	>	RAY     *rp;            /* originating ray sample */
26	>	FVECT   ux, uy;         /* tangent axis unit vectors */
27	>	int     ns;             /* number of samples per axis */
28	>	COLOR   acoef;          /* division contribution coefficient */
29	>	struct s_ambsamp {
30	>	COLOR   v;              /* hemisphere sample value */
31	>	float   p[3];           /* intersection point */
32	>	} sa[1];                /* sample array (extends struct) */
33		}  AMBHEMI;             /* ambient sample hemisphere */
34
35	<	extern double  sin(), cos(), sqrt();
35	>	#define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
36
37	+	typedef struct {
38	+	FVECT   r_i, r_i1, e_i;
39	+	double  nf, I1, I2, J2;
40	+	} FFTRI;                /* vectors and coefficients for Hessian calculation */
41
42	<	static int
43	<	ambcmp(d1, d2)                          /* decreasing order */
44	<	AMBSAMP  *d1, *d2;
42	>
43	>	static AMBHEMI *
44	>	inithemi(                       /* initialize sampling hemisphere */
45	>	COLOR   ac,
46	>	RAY     *r,
47	>	double  wt
48	>	)
49		{
50	<	if (d1->k < d2->k)
51	<	return(1);
52	<	if (d1->k > d2->k)
53	<	return(-1);
54	<	return(0);
50	>	AMBHEMI *hp;
51	>	double  d;
52	>	int     n, i;
53	>	/* set number of divisions */
54	>	if (ambacc <= FTINY &&
55	>	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
56	>	wt = d;                 /* avoid ray termination */
57	>	n = sqrt(ambdiv * wt) + 0.5;
58	>	i = 1 + 5*(ambacc > FTINY);     /* minimum number of samples */
59	>	if (n < i)
60	>	n = i;
61	>	/* allocate sampling array */
62	>	hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
63	>	sizeof(struct s_ambsamp)*(n*n - 1));
64	>	if (hp == NULL)
65	>	return(NULL);
66	>	hp->rp = r;
67	>	hp->ns = n;
68	>	/* assign coefficient */
69	>	copycolor(hp->acoef, ac);
70	>	d = 1.0/(n*n);
71	>	scalecolor(hp->acoef, d);
72	>	/* make tangent plane axes */
73	>	hp->uy[0] = 0.1 - 0.2*frandom();
74	>	hp->uy[1] = 0.1 - 0.2*frandom();
75	>	hp->uy[2] = 0.1 - 0.2*frandom();
76	>	for (i = 0; i < 3; i++)
77	>	if (r->ron[i] < 0.6 && r->ron[i] > -0.6)
78	>	break;
79	>	if (i >= 3)
80	>	error(CONSISTENCY, "bad ray direction in inithemi()");
81	>	hp->uy[i] = 1.0;
82	>	VCROSS(hp->ux, hp->uy, r->ron);
83	>	normalize(hp->ux);
84	>	VCROSS(hp->uy, r->ron, hp->ux);
85	>	/* we're ready to sample */
86	>	return(hp);
87		}
88
89
90	+	static struct s_ambsamp *
91	+	ambsample(                              /* sample an ambient direction */
92	+	AMBHEMI *hp,
93	+	int     i,
94	+	int     j
95	+	)
96	+	{
97	+	struct s_ambsamp        *ap = &ambsamp(hp,i,j);
98	+	RAY                     ar;
99	+	int                     hlist[3];
100	+	double                  spt[2], zd;
101	+	int                     ii;
102	+	/* ambient coefficient for weight */
103	+	if (ambacc > FTINY)
104	+	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
105	+	else
106	+	copycolor(ar.rcoef, hp->acoef);
107	+	if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) {
108	+	setcolor(ap->v, 0., 0., 0.);
109	+	VCOPY(ap->p, hp->rp->rop);
110	+	return(NULL);           /* no sample taken */
111	+	}
112	+	if (ambacc > FTINY) {
113	+	multcolor(ar.rcoef, hp->acoef);
114	+	scalecolor(ar.rcoef, 1./AVGREFL);
115	+	}
116	+	/* generate hemispherical sample */
117	+	SDsquare2disk(spt,      (i+.1+.8*frandom())/hp->ns,
118	+	(j+.1+.8*frandom())/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] +
122	+	spt[1]*hp->uy[ii] +
123	+	zd*hp->rp->ron[ii];
124	+	checknorm(ar.rdir);
125	+	dimlist[ndims++] = i*hp->ns + j + 90171;
126	+	rayvalue(&ar);                  /* evaluate ray */
127	+	ndims--;
128	+	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
129	+	copycolor(ap->v, ar.rcol);
130	+	if (ar.rt > 20.0*maxarad)       /* limit vertex distance */
131	+	ar.rt = 20.0*maxarad;
132	+	VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
133	+	return(ap);
134	+	}
135	+
136	+
137	+	/* Compute vectors and coefficients for Hessian/gradient calcs */
138	+	static void
139	+	comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop)
140	+	{
141	+	FVECT   v1;
142	+	double  dot_e, dot_er, dot_r, dot_r1;
143	+
144	+	VSUB(ftp->r_i, ap0, rop);
145	+	VSUB(ftp->r_i1, ap1, rop);
146	+	VSUB(ftp->e_i, ap1, ap0);
147	+	VCROSS(v1, ftp->e_i, ftp->r_i);
148	+	ftp->nf = 1.0/DOT(v1,v1);
149	+	VCROSS(v1, ftp->r_i, ftp->r_i1);
150	+	ftp->I1 = sqrt(DOT(v1,v1)*ftp->nf);
151	+	dot_e = DOT(ftp->e_i,ftp->e_i);
152	+	dot_er = DOT(ftp->e_i, ftp->r_i);
153	+	dot_r = DOT(ftp->r_i,ftp->r_i);
154	+	dot_r1 = DOT(ftp->r_i1,ftp->r_i1);
155	+	ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r +
156	+	dot_e*ftp->I1 )*0.5*ftp->nf;
157	+	ftp->J2 =  0.25*ftp->nf*( 1.0/dot_r - 1.0/dot_r1 ) -
158	+	dot_er/dot_e*ftp->I2;
159	+	}
160	+
161	+
162	+	/* Compose 3x3 matrix from two vectors */
163	+	static void
164	+	compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
165	+	{
166	+	mat[0][0] = 2.0*va[0]*vb[0];
167	+	mat[1][1] = 2.0*va[1]*vb[1];
168	+	mat[2][2] = 2.0*va[2]*vb[2];
169	+	mat[0][1] = mat[1][0] = va[0]*vb[1] + va[1]*vb[0];
170	+	mat[0][2] = mat[2][0] = va[0]*vb[2] + va[2]*vb[0];
171	+	mat[1][2] = mat[2][1] = va[1]*vb[2] + va[2]*vb[1];
172	+	}
173	+
174	+
175	+	/* Compute partial 3x3 Hessian matrix for edge */
176	+	static void
177	+	comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
178	+	{
179	+	FVECT   v1, v2;
180	+	FVECT   m1[3], m2[3], m3[3], m4[3];
181	+	double  d1, d2, d3, d4;
182	+	double  I3, J3, K3;
183	+	int     i, j;
184	+	/* compute intermediate coefficients */
185	+	d1 = 1.0/DOT(ftp->r_i,ftp->r_i);
186	+	d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
187	+	d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
188	+	d4 = DOT(ftp->e_i, ftp->r_i);
189	+	I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 +
190	+	3.0*d3*ftp->I2 );
191	+	J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3;
192	+	K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3);
193	+	/* intermediate matrices */
194	+	VCROSS(v1, nrm, ftp->e_i);
195	+	for (j = 3; j--; )
196	+	v2[j] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j];
197	+	compose_matrix(m1, v1, v2);
198	+	compose_matrix(m2, ftp->r_i, ftp->r_i);
199	+	compose_matrix(m3, ftp->e_i, ftp->e_i);
200	+	compose_matrix(m4, ftp->r_i, ftp->e_i);
201	+	VCROSS(v1, ftp->r_i, ftp->e_i);
202	+	d1 = DOT(nrm, v1);
203	+	d2 = -d1*ftp->I2;
204	+	d1 *= 2.0;
205	+	for (i = 3; i--; )              /* final matrix sum */
206	+	for (j = 3; j--; ) {
207	+	hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] +
208	+	2.0*J3*m4[i][j] );
209	+	hess[i][j] += d2*(i==j);
210	+	hess[i][j] *= -1.0/PI;
211	+	}
212	+	}
213	+
214	+
215	+	/* Reverse hessian calculation result for edge in other direction */
216	+	static void
217	+	rev_hessian(FVECT hess[3])
218	+	{
219	+	int     i;
220	+
221	+	for (i = 3; i--; ) {
222	+	hess[i][0] = -hess[i][0];
223	+	hess[i][1] = -hess[i][1];
224	+	hess[i][2] = -hess[i][2];
225	+	}
226	+	}
227	+
228	+
229	+	/* Add to radiometric Hessian from the given triangle */
230	+	static void
231	+	add2hessian(FVECT hess[3], FVECT ehess1[3],
232	+	FVECT ehess2[3], FVECT ehess3[3], COLORV v)
233	+	{
234	+	int     i, j;
235	+
236	+	for (i = 3; i--; )
237	+	for (j = 3; j--; )
238	+	hess[i][j] += v*( ehess1[i][j] + ehess2[i][j] + ehess3[i][j] );
239	+	}
240	+
241	+
242	+	/* Compute partial displacement form factor gradient for edge */
243	+	static void
244	+	comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
245	+	{
246	+	FVECT   vcp;
247	+	double  f1;
248	+	int     i;
249	+
250	+	VCROSS(vcp, ftp->r_i, ftp->r_i1);
251	+	f1 = 2.0*DOT(nrm, vcp);
252	+	VCROSS(vcp, nrm, ftp->e_i);
253	+	for (i = 3; i--; )
254	+	grad[i] = (0.5/PI)*( ftp->I1*vcp[i] +
255	+	f1*(ftp->I2*ftp->r_i[i] + ftp->J2*ftp->e_i[i]) );
256	+	}
257	+
258	+
259	+	/* Reverse gradient calculation result for edge in other direction */
260	+	static void
261	+	rev_gradient(FVECT grad)
262	+	{
263	+	grad[0] = -grad[0];
264	+	grad[1] = -grad[1];
265	+	grad[2] = -grad[2];
266	+	}
267	+
268	+
269	+	/* Add to displacement gradient from the given triangle */
270	+	static void
271	+	add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
272	+	{
273	+	int     i;
274	+
275	+	for (i = 3; i--; )
276	+	grad[i] += v*( egrad1[i] + egrad2[i] + egrad3[i] );
277	+	}
278	+
279	+
280	+	/* Return brightness of furthest ambient sample */
281	+	static COLORV
282	+	back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2,
283	+	struct s_ambsamp *ap3, FVECT orig)
284	+	{
285	+	COLORV  vback;
286	+	FVECT   vec;
287	+	double  d2, d2best;
288	+
289	+	VSUB(vec, ap1->p, orig);
290	+	d2best = DOT(vec,vec);
291	+	vback = ap1->v[CIEY];
292	+	VSUB(vec, ap2->p, orig);
293	+	d2 = DOT(vec,vec);
294	+	if (d2 > d2best) {
295	+	d2best = d2;
296	+	vback = ap2->v[CIEY];
297	+	}
298	+	VSUB(vec, ap3->p, orig);
299	+	d2 = DOT(vec,vec);
300	+	if (d2 > d2best)
301	+	return(ap3->v[CIEY]);
302	+	return(vback);
303	+	}
304	+
305	+
306	+	/* Compute anisotropic radii and eigenvector directions */
307		static int
308	<	ambnorm(d1, d2)                         /* standard order */
47	<	AMBSAMP  *d1, *d2;
308	>	eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
309		{
310	<	register int  c;
310	>	double  hess2[2][2];
311	>	FVECT   a, b;
312	>	double  evalue[2], slope1, xmag1;
313	>	int     i;
314	>	/* project Hessian to sample plane */
315	>	for (i = 3; i--; ) {
316	>	a[i] = DOT(hessian[i], uv[0]);
317	>	b[i] = DOT(hessian[i], uv[1]);
318	>	}
319	>	hess2[0][0] = DOT(uv[0], a);
320	>	hess2[0][1] = DOT(uv[0], b);
321	>	hess2[1][0] = DOT(uv[1], a);
322	>	hess2[1][1] = DOT(uv[1], b);
323	>	/* compute eigenvalues */
324	>	if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
325	>	hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 ||
326	>	(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY ||
327	>	(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY )
328	>	error(INTERNAL, "bad eigenvalue calculation");
329
330	<	if (c = d1->t - d2->t)
331	<	return(c);
332	<	return(d1->p - d2->p);
330	>	if (evalue[0] > evalue[1]) {
331	>	ra[0] = 1.0/sqrt(sqrt(evalue[0]));
332	>	ra[1] = 1.0/sqrt(sqrt(evalue[1]));
333	>	slope1 = evalue[1];
334	>	} else {
335	>	ra[0] = 1.0/sqrt(sqrt(evalue[1]));
336	>	ra[1] = 1.0/sqrt(sqrt(evalue[0]));
337	>	slope1 = evalue[0];
338	>	}
339	>	/* compute unit eigenvectors */
340	>	if (fabs(hess2[0][1]) <= FTINY)
341	>	return;                 /* uv OK as is */
342	>	slope1 = (slope1 - hess2[0][0]) / hess2[0][1];
343	>	xmag1 = sqrt(1.0/(1.0 + slope1*slope1));
344	>	for (i = 3; i--; ) {
345	>	b[i] = xmag1*uv[0][i] + slope1*xmag1*uv[1][i];
346	>	a[i] = slope1*xmag1*uv[0][i] - xmag1*uv[1][i];
347	>	}
348	>	VCOPY(uv[0], a);
349	>	VCOPY(uv[1], b);
350		}
351
352
353	<	divsample(dp, h, r)                     /* sample a division */
354	<	register AMBSAMP  *dp;
355	<	AMBHEMI  *h;
356	<	RAY  *r;
353	>	static void
354	>	ambHessian(                             /* anisotropic radii & pos. gradient */
355	>	AMBHEMI *hp,
356	>	FVECT   uv[2],                  /* returned */
357	>	float   ra[2],                  /* returned (optional) */
358	>	float   pg[2]                   /* returned (optional) */
359	>	)
360		{
361	+	static char     memerrmsg[] = "out of memory in ambHessian()";
362	+	FVECT           (*hessrow)[3] = NULL;
363	+	FVECT           *gradrow = NULL;
364	+	FVECT           hessian[3];
365	+	FVECT           gradient;
366	+	FFTRI           fftr;
367	+	int             i, j;
368	+	/* be sure to assign unit vectors */
369	+	VCOPY(uv[0], hp->ux);
370	+	VCOPY(uv[1], hp->uy);
371	+	/* clock-wise vertex traversal from sample POV */
372	+	if (ra != NULL) {               /* initialize Hessian row buffer */
373	+	hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*(hp->ns-1));
374	+	if (hessrow == NULL)
375	+	error(SYSTEM, memerrmsg);
376	+	memset(hessian, 0, sizeof(hessian));
377	+	} else if (pg == NULL)          /* bogus call? */
378	+	return;
379	+	if (pg != NULL) {               /* initialize form factor row buffer */
380	+	gradrow = (FVECT *)malloc(sizeof(FVECT)*(hp->ns-1));
381	+	if (gradrow == NULL)
382	+	error(SYSTEM, memerrmsg);
383	+	memset(gradient, 0, sizeof(gradient));
384	+	}
385	+	/* compute first row of edges */
386	+	for (j = 0; j < hp->ns-1; j++) {
387	+	comp_fftri(&fftr, ambsamp(hp,0,j).p,
388	+	ambsamp(hp,0,j+1).p, hp->rp->rop);
389	+	if (hessrow != NULL)
390	+	comp_hessian(hessrow[j], &fftr, hp->rp->ron);
391	+	if (gradrow != NULL)
392	+	comp_gradient(gradrow[j], &fftr, hp->rp->ron);
393	+	}
394	+	/* sum each row of triangles */
395	+	for (i = 0; i < hp->ns-1; i++) {
396	+	FVECT       hesscol[3];     /* compute first vertical edge */
397	+	FVECT       gradcol;
398	+	comp_fftri(&fftr, ambsamp(hp,i,0).p,
399	+	ambsamp(hp,i+1,0).p, hp->rp->rop);
400	+	if (hessrow != NULL)
401	+	comp_hessian(hesscol, &fftr, hp->rp->ron);
402	+	if (gradrow != NULL)
403	+	comp_gradient(gradcol, &fftr, hp->rp->ron);
404	+	for (j = 0; j < hp->ns-1; j++) {
405	+	FVECT   hessdia[3];     /* compute triangle contributions */
406	+	FVECT   graddia;
407	+	COLORV  backg;
408	+	backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
409	+	&ambsamp(hp,i+1,j), hp->rp->rop);
410	+	/* diagonal (inner) edge */
411	+	comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
412	+	ambsamp(hp,i+1,j).p, hp->rp->rop);
413	+	if (hessrow != NULL) {
414	+	comp_hessian(hessdia, &fftr, hp->rp->ron);
415	+	rev_hessian(hesscol);
416	+	add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
417	+	}
418	+	if (gradient != NULL) {
419	+	comp_gradient(graddia, &fftr, hp->rp->ron);
420	+	rev_gradient(gradcol);
421	+	add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
422	+	}
423	+	/* initialize edge in next row */
424	+	comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
425	+	ambsamp(hp,i+1,j).p, hp->rp->rop);
426	+	if (hessrow != NULL)
427	+	comp_hessian(hessrow[j], &fftr, hp->rp->ron);
428	+	if (gradrow != NULL)
429	+	comp_gradient(gradrow[j], &fftr, hp->rp->ron);
430	+	/* new column edge & paired triangle */
431	+	backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
432	+	&ambsamp(hp,i+1,j), hp->rp->rop);
433	+	comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
434	+	hp->rp->rop);
435	+	if (hessrow != NULL) {
436	+	comp_hessian(hesscol, &fftr, hp->rp->ron);
437	+	rev_hessian(hessdia);
438	+	add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
439	+	if (i < hp->ns-2)
440	+	rev_hessian(hessrow[j]);
441	+	}
442	+	if (gradrow != NULL) {
443	+	comp_gradient(gradcol, &fftr, hp->rp->ron);
444	+	rev_gradient(graddia);
445	+	add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
446	+	if (i < hp->ns-2)
447	+	rev_gradient(gradrow[j]);
448	+	}
449	+	}
450	+	}
451	+	/* release row buffers */
452	+	if (hessrow != NULL) free(hessrow);
453	+	if (gradrow != NULL) free(gradrow);
454	+
455	+	if (ra != NULL)                 /* extract eigenvectors & radii */
456	+	eigenvectors(uv, ra, hessian);
457	+	if (pg != NULL) {               /* project position gradient */
458	+	pg[0] = DOT(gradient, uv[0]);
459	+	pg[1] = DOT(gradient, uv[1]);
460	+	}
461	+	}
462	+
463	+
464	+	/* Compute direction gradient from a hemispherical sampling */
465	+	static void
466	+	ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
467	+	{
468	+	struct s_ambsamp        *ap;
469	+	int                     n;
470	+	FVECT                   vd;
471	+	double                  gfact;
472	+
473	+	dg[0] = dg[1] = 0;
474	+	for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
475	+	/* use vector for azimuth + 90deg */
476	+	VSUB(vd, ap->p, hp->rp->rop);
477	+	/* brightness with tangent factor */
478	+	gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd);
479	+	/* sine = proj_radius/vd_length */
480	+	dg[0] -= DOT(uv[1], vd) * gfact;
481	+	dg[1] += DOT(uv[0], vd) * gfact;
482	+	}
483	+	}
484	+
485	+
486	+	int
487	+	doambient(                              /* compute ambient component */
488	+	COLOR   rcol,                   /* input/output color */
489	+	RAY     *r,
490	+	double  wt,
491	+	FVECT   uv[2],                  /* returned (optional) */
492	+	float   ra[2],                  /* returned (optional) */
493	+	float   pg[2],                  /* returned (optional) */
494	+	float   dg[2]                   /* returned (optional) */
495	+	)
496	+	{
497	+	AMBHEMI                 *hp = inithemi(rcol, r, wt);
498	+	int                     cnt = 0;
499	+	FVECT                   my_uv[2];
500	+	double                  d, acol[3];
501	+	struct s_ambsamp        *ap;
502	+	int                     i, j;
503	+	/* check/initialize */
504	+	if (hp == NULL)
505	+	return(0);
506	+	if (uv != NULL)
507	+	memset(uv, 0, sizeof(FVECT)*2);
508	+	if (ra != NULL)
509	+	ra[0] = ra[1] = 0.0;
510	+	if (pg != NULL)
511	+	pg[0] = pg[1] = 0.0;
512	+	if (dg != NULL)
513	+	dg[0] = dg[1] = 0.0;
514	+	/* sample the hemisphere */
515	+	acol[0] = acol[1] = acol[2] = 0.0;
516	+	for (i = hp->ns; i--; )
517	+	for (j = hp->ns; j--; )
518	+	if ((ap = ambsample(hp, i, j)) != NULL) {
519	+	addcolor(acol, ap->v);
520	+	++cnt;
521	+	}
522	+	if (!cnt) {
523	+	setcolor(rcol, 0.0, 0.0, 0.0);
524	+	free(hp);
525	+	return(0);              /* no valid samples */
526	+	}
527	+	d = 1.0 / cnt;                  /* final indirect irradiance/PI */
528	+	acol[0] *= d; acol[1] *= d; acol[2] *= d;
529	+	copycolor(rcol, acol);
530	+	if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
531	+	(ra == NULL) & (pg == NULL) & (dg == NULL)) {
532	+	free(hp);
533	+	return(-1);             /* no radius or gradient calc. */
534	+	}
535	+	d = 0.01 * bright(rcol);        /* add in 1% before Hessian comp. */
536	+	if (d < FTINY) d = FTINY;
537	+	ap = hp->sa;                    /* using Y channel from here on... */
538	+	for (i = hp->ns*hp->ns; i--; ap++)
539	+	colval(ap->v,CIEY) = bright(ap->v) + d;
540	+
541	+	if (uv == NULL)                 /* make sure we have axis pointers */
542	+	uv = my_uv;
543	+	/* compute radii & pos. gradient */
544	+	ambHessian(hp, uv, ra, pg);
545	+	if (dg != NULL)                 /* compute direction gradient */
546	+	ambdirgrad(hp, uv, dg);
547	+	if (ra != NULL) {               /* scale/clamp radii */
548	+	d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt));
549	+	ra[0] *= d;
550	+	if ((ra[1] *= d) > 2.0*ra[0])
551	+	ra[1] = 2.0*ra[0];
552	+	if (ra[1] > maxarad) {
553	+	ra[1] = maxarad;
554	+	if (ra[0] > maxarad)
555	+	ra[0] = maxarad;
556	+	}
557	+	}
558	+	free(hp);                       /* clean up and return */
559	+	return(1);
560	+	}
561	+
562	+
563	+	#else /* ! NEWAMB */
564	+
565	+
566	+	void
567	+	inithemi(                       /* initialize sampling hemisphere */
568	+	AMBHEMI  *hp,
569	+	COLOR ac,
570	+	RAY  *r,
571	+	double  wt
572	+	)
573	+	{
574	+	double  d;
575	+	int  i;
576	+	/* set number of divisions */
577	+	if (ambacc <= FTINY &&
578	+	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
579	+	wt = d;                 /* avoid ray termination */
580	+	hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
581	+	i = ambacc > FTINY ? 3 : 1;     /* minimum number of samples */
582	+	if (hp->nt < i)
583	+	hp->nt = i;
584	+	hp->np = PI * hp->nt + 0.5;
585	+	/* set number of super-samples */
586	+	hp->ns = ambssamp * wt + 0.5;
587	+	/* assign coefficient */
588	+	copycolor(hp->acoef, ac);
589	+	d = 1.0/(hp->nt*hp->np);
590	+	scalecolor(hp->acoef, d);
591	+	/* make axes */
592	+	VCOPY(hp->uz, r->ron);
593	+	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
594	+	for (i = 0; i < 3; i++)
595	+	if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
596	+	break;
597	+	if (i >= 3)
598	+	error(CONSISTENCY, "bad ray direction in inithemi");
599	+	hp->uy[i] = 1.0;
600	+	fcross(hp->ux, hp->uy, hp->uz);
601	+	normalize(hp->ux);
602	+	fcross(hp->uy, hp->uz, hp->ux);
603	+	}
604	+
605	+
606	+	int
607	+	divsample(                              /* sample a division */
608	+	AMBSAMP  *dp,
609	+	AMBHEMI  *h,
610	+	RAY  *r
611	+	)
612	+	{
613		RAY  ar;
614	<	int  hlist[4];
614	>	int  hlist[3];
615	>	double  spt[2];
616		double  xd, yd, zd;
617		double  b2;
618		double  phi;
619	<	register int  i;
620	<
621	<	if (rayorigin(&ar, r, AMBIENT, 0.5) < 0)
619	>	int  i;
620	>	/* ambient coefficient for weight */
621	>	if (ambacc > FTINY)
622	>	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
623	>	else
624	>	copycolor(ar.rcoef, h->acoef);
625	>	if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
626		return(-1);
627	+	if (ambacc > FTINY) {
628	+	multcolor(ar.rcoef, h->acoef);
629	+	scalecolor(ar.rcoef, 1./AVGREFL);
630	+	}
631		hlist[0] = r->rno;
632		hlist[1] = dp->t;
633		hlist[2] = dp->p;
634	<	hlist[3] = 0;
635	<	zd = sqrt((dp->t+urand(ilhash(hlist,4)+dp->n))/h->nt);
636	<	hlist[3] = 1;
637	<	phi = 2.0*PI * (dp->p+urand(ilhash(hlist,4)+dp->n))/h->np;
638	<	xd = cos(phi) * zd;
79	<	yd = sin(phi) * zd;
634	>	multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
635	>	zd = sqrt((dp->t + spt[0])/h->nt);
636	>	phi = 2.0*PI * (dp->p + spt[1])/h->np;
637	>	xd = tcos(phi) * zd;
638	>	yd = tsin(phi) * zd;
639		zd = sqrt(1.0 - zd*zd);
640		for (i = 0; i < 3; i++)
641		ar.rdir[i] =    xd*h->ux[i] +
642		yd*h->uy[i] +
643		zd*h->uz[i];
644	+	checknorm(ar.rdir);
645		dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
646		rayvalue(&ar);
647		ndims--;
648	+	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
649		addcolor(dp->v, ar.rcol);
650	<	if (ar.rt < FHUGE)
650	>	/* use rt to improve gradient calc */
651	>	if (ar.rt > FTINY && ar.rt < FHUGE)
652		dp->r += 1.0/ar.rt;
653		/* (re)initialize error */
654		if (dp->n++) {
#	Line 99 | Line 661 | RAY  *r;
661		}
662
663
664	+	static int
665	+	ambcmp(                                 /* decreasing order */
666	+	const void *p1,
667	+	const void *p2
668	+	)
669	+	{
670	+	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
671	+	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
672	+
673	+	if (d1->k < d2->k)
674	+	return(1);
675	+	if (d1->k > d2->k)
676	+	return(-1);
677	+	return(0);
678	+	}
679	+
680	+
681	+	static int
682	+	ambnorm(                                /* standard order */
683	+	const void *p1,
684	+	const void *p2
685	+	)
686	+	{
687	+	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
688	+	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
689	+	int     c;
690	+
691	+	if ( (c = d1->t - d2->t) )
692	+	return(c);
693	+	return(d1->p - d2->p);
694	+	}
695	+
696	+
697		double
698	<	doambient(acol, r, pg, dg)              /* compute ambient component */
699	<	COLOR  acol;
700	<	RAY  *r;
701	<	FVECT  pg, dg;
698	>	doambient(                              /* compute ambient component */
699	>	COLOR  rcol,
700	>	RAY  *r,
701	>	double  wt,
702	>	FVECT  pg,
703	>	FVECT  dg
704	>	)
705		{
706	<	double  b, d;
706	>	double  b, d=0;
707		AMBHEMI  hemi;
708		AMBSAMP  *div;
709		AMBSAMP  dnew;
710	<	register AMBSAMP  *dp;
710	>	double  acol[3];
711	>	AMBSAMP  *dp;
712		double  arad;
713	<	int  ndivs, ns;
714	<	register int  i, j;
116	<	/* initialize color */
117	<	setcolor(acol, 0.0, 0.0, 0.0);
713	>	int  divcnt;
714	>	int  i, j;
715		/* initialize hemisphere */
716	<	inithemi(&hemi, r);
717	<	ndivs = hemi.nt * hemi.np;
718	<	if (ndivs == 0)
716	>	inithemi(&hemi, rcol, r, wt);
717	>	divcnt = hemi.nt * hemi.np;
718	>	/* initialize */
719	>	if (pg != NULL)
720	>	pg[0] = pg[1] = pg[2] = 0.0;
721	>	if (dg != NULL)
722	>	dg[0] = dg[1] = dg[2] = 0.0;
723	>	setcolor(rcol, 0.0, 0.0, 0.0);
724	>	if (divcnt == 0)
725		return(0.0);
726	<	/* set number of super-samples */
727	<	ns = ambssamp * r->rweight + 0.5;
728	<	if (ns > 0 || pg != NULL || dg != NULL) {
126	<	div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
726	>	/* allocate super-samples */
727	>	if (hemi.ns > 0 || pg != NULL || dg != NULL) {
728	>	div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP));
729		if (div == NULL)
730		error(SYSTEM, "out of memory in doambient");
731		} else
732		div = NULL;
733		/* sample the divisions */
734		arad = 0.0;
735	+	acol[0] = acol[1] = acol[2] = 0.0;
736		if ((dp = div) == NULL)
737		dp = &dnew;
738	+	divcnt = 0;
739		for (i = 0; i < hemi.nt; i++)
740		for (j = 0; j < hemi.np; j++) {
741		dp->t = i; dp->p = j;
742		setcolor(dp->v, 0.0, 0.0, 0.0);
743		dp->r = 0.0;
744		dp->n = 0;
745	<	if (divsample(dp, &hemi, r) < 0)
746	<	goto oopsy;
745	>	if (divsample(dp, &hemi, r) < 0) {
746	>	if (div != NULL)
747	>	dp++;
748	>	continue;
749	>	}
750	>	arad += dp->r;
751	>	divcnt++;
752		if (div != NULL)
753		dp++;
754	<	else {
754	>	else
755		addcolor(acol, dp->v);
147	–	arad += dp->r;
148	–	}
756		}
757	<	if (ns > 0) {                   /* perform super-sampling */
757	>	if (!divcnt) {
758	>	if (div != NULL)
759	>	free((void *)div);
760	>	return(0.0);            /* no samples taken */
761	>	}
762	>	if (divcnt < hemi.nt*hemi.np) {
763	>	pg = dg = NULL;         /* incomplete sampling */
764	>	hemi.ns = 0;
765	>	} else if (arad > FTINY && divcnt/arad < minarad) {
766	>	hemi.ns = 0;            /* close enough */
767	>	} else if (hemi.ns > 0) {       /* else perform super-sampling? */
768		comperrs(div, &hemi);                   /* compute errors */
769	<	qsort(div, ndivs, sizeof(AMBSAMP), ambcmp);     /* sort divs */
769	>	qsort(div, divcnt, sizeof(AMBSAMP), ambcmp);    /* sort divs */
770		/* super-sample */
771	<	for (i = ns; i > 0; i--) {
772	<	copystruct(&dnew, div);
773	<	if (divsample(&dnew, &hemi, r) < 0)
774	<	goto oopsy;
775	<	/* reinsert */
776	<	dp = div;
777	<	j = ndivs < i ? ndivs : i;
771	>	for (i = hemi.ns; i > 0; i--) {
772	>	dnew = *div;
773	>	if (divsample(&dnew, &hemi, r) < 0) {
774	>	dp++;
775	>	continue;
776	>	}
777	>	dp = div;               /* reinsert */
778	>	j = divcnt < i ? divcnt : i;
779		while (--j > 0 && dnew.k < dp[1].k) {
780	<	copystruct(dp, dp+1);
780	>	*dp = *(dp+1);
781		dp++;
782		}
783	<	copystruct(dp, &dnew);
783	>	*dp = dnew;
784		}
785		if (pg != NULL || dg != NULL)   /* restore order */
786	<	qsort(div, ndivs, sizeof(AMBSAMP), ambnorm);
786	>	qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
787		}
788		/* compute returned values */
789		if (div != NULL) {
790	<	for (i = ndivs, dp = div; i-- > 0; dp++) {
790	>	arad = 0.0;             /* note: divcnt may be < nt*np */
791	>	for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
792		arad += dp->r;
793		if (dp->n > 1) {
794		b = 1.0/dp->n;
#	Line 179 | Line 798 | FVECT  pg, dg;
798		}
799		addcolor(acol, dp->v);
800		}
801	<	if (pg != NULL)
802	<	posgradient(pg, div, &hemi);
803	<	if (dg != NULL)
804	<	dirgradient(dg, div, &hemi);
805	<	free((char *)div);
801	>	b = bright(acol);
802	>	if (b > FTINY) {
803	>	b = 1.0/b;      /* compute & normalize gradient(s) */
804	>	if (pg != NULL) {
805	>	posgradient(pg, div, &hemi);
806	>	for (i = 0; i < 3; i++)
807	>	pg[i] *= b;
808	>	}
809	>	if (dg != NULL) {
810	>	dirgradient(dg, div, &hemi);
811	>	for (i = 0; i < 3; i++)
812	>	dg[i] *= b;
813	>	}
814	>	}
815	>	free((void *)div);
816		}
817	<	b = 1.0/ndivs;
189	<	scalecolor(acol, b);
817	>	copycolor(rcol, acol);
818		if (arad <= FTINY)
191	–	arad = FHUGE;
192	–	else
193	–	arad = (ndivs+ns)/arad;
194	–	if (arad > maxarad)
819		arad = maxarad;
820	<	else if (arad < minarad)
820	>	else
821	>	arad = (divcnt+hemi.ns)/arad;
822	>	if (pg != NULL) {               /* reduce radius if gradient large */
823	>	d = DOT(pg,pg);
824	>	if (d*arad*arad > 1.0)
825	>	arad = 1.0/sqrt(d);
826	>	}
827	>	if (arad < minarad) {
828		arad = minarad;
829	<	arad /= sqrt(r->rweight);
829	>	if (pg != NULL && d*arad*arad > 1.0) {  /* cap gradient */
830	>	d = 1.0/arad/sqrt(d);
831	>	for (i = 0; i < 3; i++)
832	>	pg[i] *= d;
833	>	}
834	>	}
835	>	if ((arad /= sqrt(wt)) > maxarad)
836	>	arad = maxarad;
837		return(arad);
200	–	oopsy:
201	–	if (div != NULL)
202	–	free((char *)div);
203	–	return(0.0);
838		}
839
840
841	<	inithemi(hp, r)                 /* initialize sampling hemisphere */
842	<	register AMBHEMI  *hp;
843	<	RAY  *r;
841	>	void
842	>	comperrs(                       /* compute initial error estimates */
843	>	AMBSAMP  *da,   /* assumes standard ordering */
844	>	AMBHEMI  *hp
845	>	)
846		{
211	–	register int  i;
212	–	/* set number of divisions */
213	–	hp->nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5;
214	–	hp->np = 2 * hp->nt;
215	–	/* make axes */
216	–	VCOPY(hp->uz, r->ron);
217	–	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
218	–	for (i = 0; i < 3; i++)
219	–	if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
220	–	break;
221	–	if (i >= 3)
222	–	error(CONSISTENCY, "bad ray direction in inithemi");
223	–	hp->uy[i] = 1.0;
224	–	fcross(hp->ux, hp->uy, hp->uz);
225	–	normalize(hp->ux);
226	–	fcross(hp->uy, hp->uz, hp->ux);
227	–	}
228	–
229	–
230	–	comperrs(da, hp)                /* compute initial error estimates */
231	–	AMBSAMP  *da;           /* assumes standard ordering */
232	–	register AMBHEMI  *hp;
233	–	{
847		double  b, b2;
848		int  i, j;
849	<	register AMBSAMP  *dp;
849	>	AMBSAMP  *dp;
850		/* sum differences from neighbors */
851		dp = da;
852		for (i = 0; i < hp->nt; i++)
853		for (j = 0; j < hp->np; j++) {
854	+	#ifdef  DEBUG
855	+	if (dp->t != i || dp->p != j)
856	+	error(CONSISTENCY,
857	+	"division order in comperrs");
858	+	#endif
859		b = bright(dp[0].v);
860		if (i > 0) {            /* from above */
861		b2 = bright(dp[-hp->np].v) - b;
#	Line 272 | Line 890 | register AMBHEMI  *hp;
890		}
891
892
893	<	posgradient(gv, da, hp)                         /* compute position gradient */
894	<	FVECT  gv;
895	<	AMBSAMP  *da;                   /* assumes standard ordering */
896	<	AMBHEMI  *hp;
893	>	void
894	>	posgradient(                                    /* compute position gradient */
895	>	FVECT  gv,
896	>	AMBSAMP  *da,                   /* assumes standard ordering */
897	>	AMBHEMI  *hp
898	>	)
899		{
900	<	register int  i, j;
901	<	double  b, d;
900	>	int  i, j;
901	>	double  nextsine, lastsine, b, d;
902		double  mag0, mag1;
903		double  phi, cosp, sinp, xd, yd;
904	<	register AMBSAMP  *dp;
904	>	AMBSAMP  *dp;
905
906		xd = yd = 0.0;
907		for (j = 0; j < hp->np; j++) {
908		dp = da + j;
909		mag0 = mag1 = 0.0;
910	+	lastsine = 0.0;
911		for (i = 0; i < hp->nt; i++) {
912		#ifdef  DEBUG
913		if (dp->t != i || dp->p != j)
#	Line 297 | Line 918 | AMBHEMI  *hp;
918		if (i > 0) {
919		d = dp[-hp->np].r;
920		if (dp[0].r > d) d = dp[0].r;
921	<	d *= 1.0 - sqrt((double)i/hp->nt);
921	>	/* sin(t)*cos(t)^2 */
922	>	d *= lastsine * (1.0 - (double)i/hp->nt);
923		mag0 += d*(b - bright(dp[-hp->np].v));
924		}
925	+	nextsine = sqrt((double)(i+1)/hp->nt);
926		if (j > 0) {
927		d = dp[-1].r;
928		if (dp[0].r > d) d = dp[0].r;
929	<	mag1 += d*(b - bright(dp[-1].v));
929	>	mag1 += d * (nextsine - lastsine) *
930	>	(b - bright(dp[-1].v));
931		} else {
932		d = dp[hp->np-1].r;
933		if (dp[0].r > d) d = dp[0].r;
934	<	mag1 += d*(b - bright(dp[hp->np-1].v));
934	>	mag1 += d * (nextsine - lastsine) *
935	>	(b - bright(dp[hp->np-1].v));
936		}
937		dp += hp->np;
938	+	lastsine = nextsine;
939		}
940	<	if (hp->nt > 1) {
315	<	mag0 /= (double)hp->np;
316	<	mag1 /= (double)hp->nt;
317	<	}
940	>	mag0 *= 2.0*PI / hp->np;
941		phi = 2.0*PI * (double)j/hp->np;
942	<	cosp = cos(phi); sinp = sin(phi);
942	>	cosp = tcos(phi); sinp = tsin(phi);
943		xd += mag0*cosp - mag1*sinp;
944		yd += mag0*sinp + mag1*cosp;
945		}
946		for (i = 0; i < 3; i++)
947	<	gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
947	>	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI;
948		}
949
950
951	<	dirgradient(gv, da, hp)                         /* compute direction gradient */
952	<	FVECT  gv;
953	<	AMBSAMP  *da;                   /* assumes standard ordering */
954	<	AMBHEMI  *hp;
951	>	void
952	>	dirgradient(                                    /* compute direction gradient */
953	>	FVECT  gv,
954	>	AMBSAMP  *da,                   /* assumes standard ordering */
955	>	AMBHEMI  *hp
956	>	)
957		{
958	<	register int  i, j;
958	>	int  i, j;
959		double  mag;
960		double  phi, xd, yd;
961	<	register AMBSAMP  *dp;
961	>	AMBSAMP  *dp;
962
963		xd = yd = 0.0;
964		for (j = 0; j < hp->np; j++) {
#	Line 345 | Line 970 | AMBHEMI  *hp;
970		error(CONSISTENCY,
971		"division order in dirgradient");
972		#endif
973	<	mag += sqrt((i+.5)/hp->nt)*bright(dp->v);
973	>	/* tan(t) */
974	>	mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
975		dp += hp->np;
976		}
977		phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
978	<	xd += mag * cos(phi);
979	<	yd += mag * sin(phi);
978	>	xd += mag * tcos(phi);
979	>	yd += mag * tsin(phi);
980		}
981		for (i = 0; i < 3; i++)
982	<	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])/(hp->nt*hp->np);
982	>	gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
983		}
984	+
985	+	#endif  /* ! NEWAMB */

Diff Legend

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