ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambcomp.c
Revision: 2.37
Committed: Sat Apr 26 05:09:54 2014 UTC (10 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.36: +7 -6 lines
Log Message: 
Went back to old axis initialization method

File Contents

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