ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambcomp.c
(Generate patch)

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.34 by greg, Thu Apr 24 23:15:10 2014 UTC vs.
Revision 2.46 by greg, Fri May 2 21:58:50 2014 UTC

# Line 8 | Line 8 | static const char      RCSid[] = "$Id$";
8   *      for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz
9   *      from ACM SIGGRAPH Asia 2012 conference proceedings.
10   *
11 + *  Added book-keeping optimization to avoid calculations that would
12 + *      cancel due to traversal both directions on edges that are adjacent
13 + *      to same-valued triangles.  This cuts about half of Hessian math.
14 + *
15   *  Declarations of external symbols in ambient.h
16   */
17  
# Line 21 | Line 25 | static const char      RCSid[] = "$Id$";
25  
26   extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
27  
28 +                                /* vertex direction bit positions */
29 + #define VDB_xy  0
30 + #define VDB_y   01
31 + #define VDB_x   02
32 + #define VDB_Xy  03
33 + #define VDB_xY  04
34 + #define VDB_X   05
35 + #define VDB_Y   06
36 + #define VDB_XY  07
37 +                                /* get opposite vertex direction bit */
38 + #define VDB_OPP(f)      (~(f) & 07)
39 +                                /* adjacent triangle vertex flags */
40 + static const int  adjacent_trifl[8] = {
41 +                        0,                      /* forbidden diagonal */
42 +                        1<<VDB_x|1<<VDB_y|1<<VDB_Xy,
43 +                        1<<VDB_y|1<<VDB_x|1<<VDB_xY,
44 +                        1<<VDB_y|1<<VDB_Xy|1<<VDB_X,
45 +                        1<<VDB_x|1<<VDB_xY|1<<VDB_Y,
46 +                        1<<VDB_Xy|1<<VDB_X|1<<VDB_Y,
47 +                        1<<VDB_xY|1<<VDB_Y|1<<VDB_X,
48 +                        0,                      /* forbidden diagonal */
49 +                };
50 +
51   typedef struct {
52 +        COLOR   v;              /* hemisphere sample value */
53 +        FVECT   p;              /* intersection point */
54 + } AMBSAMP;              /* sample value */
55 +
56 + typedef struct {
57          RAY     *rp;            /* originating ray sample */
58          FVECT   ux, uy;         /* tangent axis unit vectors */
59          int     ns;             /* number of samples per axis */
60          COLOR   acoef;          /* division contribution coefficient */
61 <        struct s_ambsamp {
30 <                COLOR   v;              /* hemisphere sample value */
31 <                FVECT   p;              /* intersection point */
32 <        } sa[1];                /* sample array (extends struct) */
61 >        AMBSAMP sa[1];          /* sample array (extends struct) */
62   }  AMBHEMI;             /* ambient sample hemisphere */
63  
64 < #define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
64 > #define ambndx(h,i,j)   ((i)*(h)->ns + (j))
65 > #define ambsam(h,i,j)   (h)->sa[ambndx(h,i,j)]
66  
67   typedef struct {
68 <        FVECT   r_i, r_i1, e_i, rI2_eJ2;
69 <        double  nf, I1, I2;
68 >        FVECT   r_i, r_i1, e_i, rcp, rI2_eJ2;
69 >        double  I1, I2;
70 >        int     valid;
71   } FFTRI;                /* vectors and coefficients for Hessian calculation */
72  
73  
74 + /* Get index for adjacent vertex */
75 + static int
76 + adjacent_verti(AMBHEMI *hp, int i, int j, int dbit)
77 + {
78 +        int     i0 = i*hp->ns + j;
79 +
80 +        switch (dbit) {
81 +        case VDB_y:     return(i0 - hp->ns);
82 +        case VDB_x:     return(i0 - 1);
83 +        case VDB_Xy:    return(i0 - hp->ns + 1);
84 +        case VDB_xY:    return(i0 + hp->ns - 1);
85 +        case VDB_X:     return(i0 + 1);
86 +        case VDB_Y:     return(i0 + hp->ns);
87 +                                /* the following should never occur */
88 +        case VDB_xy:    return(i0 - hp->ns - 1);
89 +        case VDB_XY:    return(i0 + hp->ns + 1);
90 +        }
91 +        return(-1);
92 + }
93 +
94 +
95 + /* Get vertex direction bit for the opposite edge to complete triangle */
96 + static int
97 + vdb_edge(int db1, int db2)
98 + {
99 +        switch (db1) {
100 +        case VDB_x:     return(db2==VDB_y ? VDB_Xy : VDB_Y);
101 +        case VDB_y:     return(db2==VDB_x ? VDB_xY : VDB_X);
102 +        case VDB_X:     return(db2==VDB_Xy ? VDB_y : VDB_xY);
103 +        case VDB_Y:     return(db2==VDB_xY ? VDB_x : VDB_Xy);
104 +        case VDB_xY:    return(db2==VDB_x ? VDB_y : VDB_X);
105 +        case VDB_Xy:    return(db2==VDB_y ? VDB_x : VDB_Y);
106 +        }
107 +        error(INTERNAL, "forbidden diagonal in vdb_edge()");
108 +        return(-1);
109 + }
110 +
111 +
112   static AMBHEMI *
113   inithemi(                       /* initialize sampling hemisphere */
114          COLOR   ac,
# Line 59 | Line 128 | inithemi(                      /* initialize sampling hemisphere */
128          if (n < i)
129                  n = i;
130                                          /* allocate sampling array */
131 <        hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
63 <                                sizeof(struct s_ambsamp)*(n*n - 1));
131 >        hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1));
132          if (hp == NULL)
133                  return(NULL);
134          hp->rp = r;
# Line 70 | Line 138 | inithemi(                      /* initialize sampling hemisphere */
138          d = 1.0/(n*n);
139          scalecolor(hp->acoef, d);
140                                          /* make tangent plane axes */
141 <        hp->uy[0] = 0.1 - 0.2*frandom();
142 <        hp->uy[1] = 0.1 - 0.2*frandom();
143 <        hp->uy[2] = 0.1 - 0.2*frandom();
144 <        for (i = 0; i < 3; i++)
145 <                if (r->ron[i] < 0.6 && r->ron[i] > -0.6)
141 >        hp->uy[0] = 0.5 - frandom();
142 >        hp->uy[1] = 0.5 - frandom();
143 >        hp->uy[2] = 0.5 - frandom();
144 >        for (i = 3; i--; )
145 >                if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6))
146                          break;
147 <        if (i >= 3)
148 <                error(CONSISTENCY, "bad ray direction in inithemi()");
147 >        if (i < 0)
148 >                error(CONSISTENCY, "bad ray direction in inithemi");
149          hp->uy[i] = 1.0;
150          VCROSS(hp->ux, hp->uy, r->ron);
151          normalize(hp->ux);
# Line 87 | Line 155 | inithemi(                      /* initialize sampling hemisphere */
155   }
156  
157  
158 < static struct s_ambsamp *
159 < ambsample(                              /* sample an ambient direction */
160 <        AMBHEMI *hp,
93 <        int     i,
94 <        int     j
95 < )
158 > /* Sample ambient division and apply weighting coefficient */
159 > static int
160 > getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
161   {
162 <        struct s_ambsamp        *ap = &ambsamp(hp,i,j);
163 <        RAY                     ar;
99 <        double                  spt[2], zd;
100 <        int                     ii;
162 >        int     hlist[3], ii;
163 >        double  spt[2], zd;
164                                          /* ambient coefficient for weight */
165          if (ambacc > FTINY)
166 <                setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
166 >                setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL);
167          else
168 <                copycolor(ar.rcoef, hp->acoef);
169 <        if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0)
170 <                goto badsample;
168 >                copycolor(arp->rcoef, hp->acoef);
169 >        if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0)
170 >                return(0);
171          if (ambacc > FTINY) {
172 <                multcolor(ar.rcoef, hp->acoef);
173 <                scalecolor(ar.rcoef, 1./AVGREFL);
172 >                multcolor(arp->rcoef, hp->acoef);
173 >                scalecolor(arp->rcoef, 1./AVGREFL);
174          }
175 <                                        /* generate hemispherical sample */
176 <        SDsquare2disk(spt,      (i+.1+.8*frandom())/hp->ns,
177 <                                (j+.1+.8*frandom())/hp->ns );
175 >        hlist[0] = hp->rp->rno;
176 >        hlist[1] = j;
177 >        hlist[2] = i;
178 >        multisamp(spt, 2, urand(ilhash(hlist,3)+n));
179 >        if (!n) {                       /* avoid border samples for n==0 */
180 >                if ((spt[0] < 0.1) | (spt[0] >= 0.9))
181 >                        spt[0] = 0.1 + 0.8*frandom();
182 >                if ((spt[1] < 0.1) | (spt[1] >= 0.9))
183 >                        spt[1] = 0.1 + 0.8*frandom();
184 >        }
185 >        SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
186          zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
187          for (ii = 3; ii--; )
188 <                ar.rdir[ii] =   spt[0]*hp->ux[ii] +
188 >                arp->rdir[ii] = spt[0]*hp->ux[ii] +
189                                  spt[1]*hp->uy[ii] +
190                                  zd*hp->rp->ron[ii];
191 <        checknorm(ar.rdir);
192 <        dimlist[ndims++] = i*hp->ns + j + 90171;
193 <        rayvalue(&ar);                  /* evaluate ray */
194 <        ndims--;
191 >        checknorm(arp->rdir);
192 >        dimlist[ndims++] = ambndx(hp,i,j) + 90171;
193 >        rayvalue(arp);                  /* evaluate ray */
194 >        ndims--;                        /* apply coefficient */
195 >        multcolor(arp->rcol, arp->rcoef);
196 >        return(1);
197 > }
198 >
199 >
200 > static AMBSAMP *
201 > ambsample(                              /* initial ambient division sample */
202 >        AMBHEMI *hp,
203 >        int     i,
204 >        int     j
205 > )
206 > {
207 >        AMBSAMP *ap = &ambsam(hp,i,j);
208 >        RAY     ar;
209 >                                        /* generate hemispherical sample */
210 >        if (!getambsamp(&ar, hp, i, j, 0))
211 >                goto badsample;
212                                          /* limit vertex distance */
213          if (ar.rt > 10.0*thescene.cusize)
214                  ar.rt = 10.0*thescene.cusize;
215          else if (ar.rt <= FTINY)        /* should never happen! */
216                  goto badsample;
217          VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
130 –        multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
218          copycolor(ap->v, ar.rcol);
219          return(ap);
220   badsample:
# Line 137 | Line 224 | badsample:
224   }
225  
226  
227 + /* Estimate errors based on ambient division differences */
228 + static float *
229 + getambdiffs(AMBHEMI *hp)
230 + {
231 +        float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
232 +        float   *ep;
233 +        AMBSAMP *ap;
234 +        double  b, d2;
235 +        int     i, j;
236 +
237 +        if (earr == NULL)               /* out of memory? */
238 +                return(NULL);
239 +                                        /* compute squared neighbor diffs */
240 +        for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
241 +            for (j = 0; j < hp->ns; j++, ap++, ep++) {
242 +                b = bright(ap[0].v);
243 +                if (i) {                /* from above */
244 +                        d2 = b - bright(ap[-hp->ns].v);
245 +                        d2 *= d2;
246 +                        ep[0] += d2;
247 +                        ep[-hp->ns] += d2;
248 +                }
249 +                if (j) {                /* from behind */
250 +                        d2 = b - bright(ap[-1].v);
251 +                        d2 *= d2;
252 +                        ep[0] += d2;
253 +                        ep[-1] += d2;
254 +                }
255 +            }
256 +                                        /* correct for number of neighbors */
257 +        earr[0] *= 2.f;
258 +        earr[hp->ns-1] *= 2.f;
259 +        earr[(hp->ns-1)*hp->ns] *= 2.f;
260 +        earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f;
261 +        for (i = 1; i < hp->ns-1; i++) {
262 +                earr[i*hp->ns] *= 4./3.;
263 +                earr[i*hp->ns + hp->ns-1] *= 4./3.;
264 +        }
265 +        for (j = 1; j < hp->ns-1; j++) {
266 +                earr[j] *= 4./3.;
267 +                earr[(hp->ns-1)*hp->ns + j] *= 4./3.;
268 +        }
269 +        return(earr);
270 + }
271 +
272 +
273 + /* Perform super-sampling on hemisphere (introduces bias) */
274 + static void
275 + ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
276 + {
277 +        float   *earr = getambdiffs(hp);
278 +        double  e2sum = 0;
279 +        AMBSAMP *ap;
280 +        RAY     ar;
281 +        COLOR   asum;
282 +        float   *ep;
283 +        int     i, j, n;
284 +
285 +        if (earr == NULL)               /* just skip calc. if no memory */
286 +                return;
287 +                                        /* add up estimated variances */
288 +        for (ep = earr + hp->ns*hp->ns; ep-- > earr; )
289 +                e2sum += *ep;
290 +        ep = earr;                      /* perform super-sampling */
291 +        for (ap = hp->sa, i = 0; i < hp->ns; i++)
292 +            for (j = 0; j < hp->ns; j++, ap++) {
293 +                int     nss = *ep/e2sum*cnt + frandom();
294 +                setcolor(asum, 0., 0., 0.);
295 +                for (n = 1; n <= nss; n++) {
296 +                        if (!getambsamp(&ar, hp, i, j, n)) {
297 +                                nss = n-1;
298 +                                break;
299 +                        }
300 +                        addcolor(asum, ar.rcol);
301 +                }
302 +                if (nss) {              /* update returned ambient value */
303 +                        const double    ssf = 1./(nss + 1);
304 +                        for (n = 3; n--; )
305 +                                acol[n] += ssf*colval(asum,n) +
306 +                                                (ssf - 1.)*colval(ap->v,n);
307 +                }
308 +                e2sum -= *ep++;         /* update remainders */
309 +                cnt -= nss;
310 +        }
311 +        free(earr);
312 + }
313 +
314 +
315 + /* Compute vertex flags, indicating farthest in each direction */
316 + static uby8 *
317 + vertex_flags(AMBHEMI *hp)
318 + {
319 +        uby8    *vflags = (uby8 *)calloc(hp->ns*hp->ns, sizeof(uby8));
320 +        double  *dist2a = (double *)malloc(sizeof(double)*hp->ns);
321 +        uby8    *vf;
322 +        int     i, j;
323 +
324 +        if ((vflags == NULL) | (dist2a == NULL))
325 +                error(SYSTEM, "out of memory in vertex_flags()");
326 +        vf = vflags;            /* compute distances along first row */
327 +        for (j = 0; j < hp->ns; j++) {
328 +                dist2a[j] = dist2(ambsam(hp,0,j).p, hp->rp->rop);
329 +                ++vf;
330 +                if (!j) continue;
331 +                if (dist2a[j] >= dist2a[j-1])
332 +                        vf[0] |= 1<<VDB_x;
333 +                else
334 +                        vf[-1] |= 1<<VDB_X;
335 +        }
336 +                                /* flag subsequent rows */
337 +        for (i = 1; i < hp->ns; i++) {
338 +            double      d2n = dist2(ambsam(hp,i,0).p, hp->rp->rop);
339 +            for (j = 0; j < hp->ns-1; j++) {
340 +                double  d2 = d2n;
341 +                if (d2 >= dist2a[j])    /* row before */
342 +                        vf[0] |= 1<<VDB_y;
343 +                else
344 +                        vf[-hp->ns] |= 1<<VDB_Y;
345 +                dist2a[j] = d2n;
346 +                if (d2 >= dist2a[j+1])  /* diagonal we care about */
347 +                        vf[0] |= 1<<VDB_Xy;
348 +                else
349 +                        vf[1-hp->ns] |= 1<<VDB_xY;
350 +                d2n = dist2(ambsam(hp,i,j+1).p, hp->rp->rop);
351 +                if (d2 >= d2n)          /* column after */
352 +                        vf[0] |= 1<<VDB_X;
353 +                else
354 +                        vf[1] |= 1<<VDB_x;
355 +                ++vf;
356 +            }
357 +            if (d2n >= dist2a[j])       /* final column edge */
358 +                vf[0] |= 1<<VDB_y;
359 +            else
360 +                vf[-hp->ns] |= 1<<VDB_Y;
361 +            dist2a[j] = d2n;
362 +            ++vf;
363 +        }
364 +        free(dist2a);
365 +        return(vflags);
366 + }
367 +
368 +
369 + /* Return brightness of farthest ambient sample */
370 + static double
371 + back_ambval(AMBHEMI *hp, int i, int j, int dbit1, int dbit2, const uby8 *vflags)
372 + {
373 +        const int       v0 = ambndx(hp,i,j);
374 +        const int       tflags = (1<<dbit1 | 1<<dbit2);
375 +        int             v1, v2;
376 +
377 +        if ((vflags[v0] & tflags) == tflags)    /* is v0 the farthest? */
378 +                return(colval(hp->sa[v0].v,CIEY));
379 +        v1 = adjacent_verti(hp, i, j, dbit1);
380 +        if (vflags[v0] & 1<<dbit2)              /* v1 farthest if v0>v2 */
381 +                return(colval(hp->sa[v1].v,CIEY));
382 +        v2 = adjacent_verti(hp, i, j, dbit2);
383 +        if (vflags[v0] & 1<<dbit1)              /* v2 farthest if v0>v1 */
384 +                return(colval(hp->sa[v2].v,CIEY));
385 +                                                /* else check if v1>v2 */
386 +        if (vflags[v1] & 1<<vdb_edge(dbit1,dbit2))
387 +                return(colval(hp->sa[v1].v,CIEY));
388 +        return(colval(hp->sa[v2].v,CIEY));
389 + }
390 +
391 +
392   /* Compute vectors and coefficients for Hessian/gradient calcs */
393   static void
394 < comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop)
394 > comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags)
395   {
396 <        FVECT   vcp;
397 <        double  dot_e, dot_er, rdot_r, rdot_r1, J2;
398 <        int     i;
396 >        const int       i0 = ambndx(hp,i,j);
397 >        double          rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
398 >        int             i1, ii;
399  
400 <        VSUB(ftp->r_i, ap0, rop);
401 <        VSUB(ftp->r_i1, ap1, rop);
402 <        VSUB(ftp->e_i, ap1, ap0);
403 <        VCROSS(vcp, ftp->e_i, ftp->r_i);
404 <        ftp->nf = 1.0/DOT(vcp,vcp);
400 >        ftp->valid = 0;                 /* check if we can skip this edge */
401 >        ii = adjacent_trifl[dbit];
402 >        if ((vflags[i0] & ii) == ii)    /* cancels if vertex used as value */
403 >                return;
404 >        i1 = adjacent_verti(hp, i, j, dbit);
405 >        ii = adjacent_trifl[VDB_OPP(dbit)];
406 >        if ((vflags[i1] & ii) == ii)    /* on either end (for both triangles) */
407 >                return;
408 >                                        /* else go ahead with calculation */
409 >        VSUB(ftp->r_i, hp->sa[i0].p, hp->rp->rop);
410 >        VSUB(ftp->r_i1, hp->sa[i1].p, hp->rp->rop);
411 >        VSUB(ftp->e_i, hp->sa[i1].p, hp->sa[i0].p);
412 >        VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1);
413 >        rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp);
414          dot_e = DOT(ftp->e_i,ftp->e_i);
415          dot_er = DOT(ftp->e_i, ftp->r_i);
416          rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i);
417          rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
418          ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) *
419 <                        sqrt( ftp->nf );
419 >                        sqrt( rdot_cp );
420          ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
421 <                        dot_e*ftp->I1 )*0.5*ftp->nf;
421 >                        dot_e*ftp->I1 )*0.5*rdot_cp;
422          J2 =  ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e;
423 <        for (i = 3; i--; )
424 <                ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i];
423 >        for (ii = 3; ii--; )
424 >                ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii];
425 >        ftp->valid++;
426   }
427  
428  
# Line 181 | Line 443 | compose_matrix(FVECT mat[3], FVECT va, FVECT vb)
443   static void
444   comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
445   {
446 <        FVECT   vcp;
446 >        FVECT   ncp;
447          FVECT   m1[3], m2[3], m3[3], m4[3];
448          double  d1, d2, d3, d4;
449          double  I3, J3, K3;
450          int     i, j;
451 +
452 +        if (!ftp->valid) {              /* preemptive test */
453 +                memset(hess, 0, sizeof(FVECT)*3);
454 +                return;
455 +        }
456                                          /* compute intermediate coefficients */
457          d1 = 1.0/DOT(ftp->r_i,ftp->r_i);
458          d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1);
459          d3 = 1.0/DOT(ftp->e_i,ftp->e_i);
460          d4 = DOT(ftp->e_i, ftp->r_i);
461 <        I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 +
462 <                                3.0/d3*ftp->I2 );
461 >        I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 )
462 >                        / ( 4.0*DOT(ftp->rcp,ftp->rcp) );
463          J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3;
464          K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3);
465                                          /* intermediate matrices */
466 <        VCROSS(vcp, nrm, ftp->e_i);
467 <        compose_matrix(m1, vcp, ftp->rI2_eJ2);
466 >        VCROSS(ncp, nrm, ftp->e_i);
467 >        compose_matrix(m1, ncp, ftp->rI2_eJ2);
468          compose_matrix(m2, ftp->r_i, ftp->r_i);
469          compose_matrix(m3, ftp->e_i, ftp->e_i);
470          compose_matrix(m4, ftp->r_i, ftp->e_i);
471 <        VCROSS(vcp, ftp->r_i, ftp->e_i);
205 <        d1 = DOT(nrm, vcp);
471 >        d1 = DOT(nrm, ftp->rcp);
472          d2 = -d1*ftp->I2;
473          d1 *= 2.0;
474          for (i = 3; i--; )              /* final matrix sum */
# Line 210 | Line 476 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
476                  hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] +
477                                                  2.0*J3*m4[i][j] );
478                  hess[i][j] += d2*(i==j);
479 <                hess[i][j] *= 1.0/PI;
479 >                hess[i][j] *= -1.0/PI;
480              }
481   }
482  
# Line 232 | Line 498 | rev_hessian(FVECT hess[3])
498   /* Add to radiometric Hessian from the given triangle */
499   static void
500   add2hessian(FVECT hess[3], FVECT ehess1[3],
501 <                FVECT ehess2[3], FVECT ehess3[3], COLORV v)
501 >                FVECT ehess2[3], FVECT ehess3[3], double v)
502   {
503          int     i, j;
504  
# Line 246 | Line 512 | add2hessian(FVECT hess[3], FVECT ehess1[3],
512   static void
513   comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
514   {
515 <        FVECT   vcp;
515 >        FVECT   ncp;
516          double  f1;
517          int     i;
518  
519 <        VCROSS(vcp, ftp->r_i, ftp->r_i1);
520 <        f1 = 2.0*DOT(nrm, vcp);
521 <        VCROSS(vcp, nrm, ftp->e_i);
519 >        if (!ftp->valid) {              /* preemptive test */
520 >                memset(grad, 0, sizeof(FVECT));
521 >                return;
522 >        }
523 >        f1 = 2.0*DOT(nrm, ftp->rcp);
524 >        VCROSS(ncp, nrm, ftp->e_i);
525          for (i = 3; i--; )
526 <                grad[i] = (-0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] );
526 >                grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] );
527   }
528  
529  
# Line 270 | Line 539 | rev_gradient(FVECT grad)
539  
540   /* Add to displacement gradient from the given triangle */
541   static void
542 < add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v)
542 > add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v)
543   {
544          int     i;
545  
# Line 279 | Line 548 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
548   }
549  
550  
282 – /* Return brightness of furthest ambient sample */
283 – static COLORV
284 – back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2,
285 –                struct s_ambsamp *ap3, FVECT orig)
286 – {
287 –        COLORV  vback;
288 –        FVECT   vec;
289 –        double  d2, d2best;
290 –
291 –        VSUB(vec, ap1->p, orig);
292 –        d2best = DOT(vec,vec);
293 –        vback = colval(ap1->v,CIEY);
294 –        VSUB(vec, ap2->p, orig);
295 –        d2 = DOT(vec,vec);
296 –        if (d2 > d2best) {
297 –                d2best = d2;
298 –                vback = colval(ap2->v,CIEY);
299 –        }
300 –        VSUB(vec, ap3->p, orig);
301 –        d2 = DOT(vec,vec);
302 –        if (d2 > d2best)
303 –                return(colval(ap3->v,CIEY));
304 –        return(vback);
305 – }
306 –
307 –
551   /* Compute anisotropic radii and eigenvector directions */
552   static int
553   eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
# Line 322 | Line 565 | eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3
565          hess2[0][1] = DOT(uv[0], b);
566          hess2[1][0] = DOT(uv[1], a);
567          hess2[1][1] = DOT(uv[1], b);
568 <                                        /* compute eigenvalues */
569 <        if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
570 <                        hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 ||
571 <                        (evalue[0] = fabs(evalue[0])) <= FTINY*FTINY ||
572 <                        (evalue[1] = fabs(evalue[1])) <= FTINY*FTINY )
568 >                                        /* compute eigenvalue(s) */
569 >        i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1],
570 >                        hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]);
571 >        if (i == 1)                     /* double-root (circle) */
572 >                evalue[1] = evalue[0];
573 >        if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) |
574 >                        ((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) )
575                  error(INTERNAL, "bad eigenvalue calculation");
576  
577          if (evalue[0] > evalue[1]) {
# Line 363 | Line 608 | ambHessian(                            /* anisotropic radii & pos. gradient */
608          static char     memerrmsg[] = "out of memory in ambHessian()";
609          FVECT           (*hessrow)[3] = NULL;
610          FVECT           *gradrow = NULL;
611 +        uby8            *vflags;
612          FVECT           hessian[3];
613          FVECT           gradient;
614          FFTRI           fftr;
# Line 384 | Line 630 | ambHessian(                            /* anisotropic radii & pos. gradient */
630                          error(SYSTEM, memerrmsg);
631                  memset(gradient, 0, sizeof(gradient));
632          }
633 +                                        /* get vertex position flags */
634 +        vflags = vertex_flags(hp);
635                                          /* compute first row of edges */
636          for (j = 0; j < hp->ns-1; j++) {
637 <                comp_fftri(&fftr, ambsamp(hp,0,j).p,
390 <                                ambsamp(hp,0,j+1).p, hp->rp->rop);
637 >                comp_fftri(&fftr, hp, 0, j, VDB_X, vflags);
638                  if (hessrow != NULL)
639                          comp_hessian(hessrow[j], &fftr, hp->rp->ron);
640                  if (gradrow != NULL)
# Line 397 | Line 644 | ambHessian(                            /* anisotropic radii & pos. gradient */
644          for (i = 0; i < hp->ns-1; i++) {
645              FVECT       hesscol[3];     /* compute first vertical edge */
646              FVECT       gradcol;
647 <            comp_fftri(&fftr, ambsamp(hp,i,0).p,
401 <                        ambsamp(hp,i+1,0).p, hp->rp->rop);
647 >            comp_fftri(&fftr, hp, i, 0, VDB_Y, vflags);
648              if (hessrow != NULL)
649                  comp_hessian(hesscol, &fftr, hp->rp->ron);
650              if (gradrow != NULL)
# Line 406 | Line 652 | ambHessian(                            /* anisotropic radii & pos. gradient */
652              for (j = 0; j < hp->ns-1; j++) {
653                  FVECT   hessdia[3];     /* compute triangle contributions */
654                  FVECT   graddia;
655 <                COLORV  backg;
656 <                backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1),
411 <                                        &ambsamp(hp,i+1,j), hp->rp->rop);
655 >                double  backg;
656 >                backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags);
657                                          /* diagonal (inner) edge */
658 <                comp_fftri(&fftr, ambsamp(hp,i,j+1).p,
414 <                                ambsamp(hp,i+1,j).p, hp->rp->rop);
658 >                comp_fftri(&fftr, hp, i, j+1, VDB_xY, vflags);
659                  if (hessrow != NULL) {
660                      comp_hessian(hessdia, &fftr, hp->rp->ron);
661                      rev_hessian(hesscol);
662                      add2hessian(hessian, hessrow[j], hessdia, hesscol, backg);
663                  }
664 <                if (gradient != NULL) {
664 >                if (gradrow != NULL) {
665                      comp_gradient(graddia, &fftr, hp->rp->ron);
666                      rev_gradient(gradcol);
667                      add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
668                  }
669                                          /* initialize edge in next row */
670 <                comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p,
427 <                                ambsamp(hp,i+1,j).p, hp->rp->rop);
670 >                comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags);
671                  if (hessrow != NULL)
672                      comp_hessian(hessrow[j], &fftr, hp->rp->ron);
673                  if (gradrow != NULL)
674                      comp_gradient(gradrow[j], &fftr, hp->rp->ron);
675                                          /* new column edge & paired triangle */
676 <                backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1),
677 <                                        &ambsamp(hp,i+1,j), hp->rp->rop);
435 <                comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p,
436 <                                hp->rp->rop);
676 >                backg = back_ambval(hp, i+1, j+1, VDB_x, VDB_y, vflags);
677 >                comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags);
678                  if (hessrow != NULL) {
679                      comp_hessian(hesscol, &fftr, hp->rp->ron);
680                      rev_hessian(hessdia);
# Line 453 | Line 694 | ambHessian(                            /* anisotropic radii & pos. gradient */
694                                          /* release row buffers */
695          if (hessrow != NULL) free(hessrow);
696          if (gradrow != NULL) free(gradrow);
697 +        free(vflags);
698          
699          if (ra != NULL)                 /* extract eigenvectors & radii */
700                  eigenvectors(uv, ra, hessian);
# Line 467 | Line 709 | ambHessian(                            /* anisotropic radii & pos. gradient */
709   static void
710   ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
711   {
712 <        struct s_ambsamp        *ap;
713 <        double                  dgsum[2];
714 <        int                     n;
715 <        FVECT                   vd;
716 <        double                  gfact;
712 >        AMBSAMP *ap;
713 >        double  dgsum[2];
714 >        int     n;
715 >        FVECT   vd;
716 >        double  gfact;
717  
718          dgsum[0] = dgsum[1] = 0.0;      /* sum values times -tan(theta) */
719          for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) {
# Line 479 | Line 721 | ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2])
721                  VSUB(vd, ap->p, hp->rp->rop);
722                                          /* brightness over cosine factor */
723                  gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd);
724 <                                        /* -sine = -proj_radius/vd_length */
725 <                dgsum[0] += DOT(uv[1], vd) * gfact;
726 <                dgsum[1] -= DOT(uv[0], vd) * gfact;
724 >                                        /* sine = proj_radius/vd_length */
725 >                dgsum[0] -= DOT(uv[1], vd) * gfact;
726 >                dgsum[1] += DOT(uv[0], vd) * gfact;
727          }
728          dg[0] = dgsum[0] / (hp->ns*hp->ns);
729          dg[1] = dgsum[1] / (hp->ns*hp->ns);
# Line 499 | Line 741 | doambient(                             /* compute ambient component */
741          float   dg[2]                   /* returned (optional) */
742   )
743   {
744 <        AMBHEMI                 *hp = inithemi(rcol, r, wt);
745 <        int                     cnt = 0;
746 <        FVECT                   my_uv[2];
747 <        double                  d, acol[3];
748 <        struct s_ambsamp        *ap;
749 <        int                     i, j;
744 >        AMBHEMI *hp = inithemi(rcol, r, wt);
745 >        int     cnt;
746 >        FVECT   my_uv[2];
747 >        double  d, K, acol[3];
748 >        AMBSAMP *ap;
749 >        int     i, j;
750                                          /* check/initialize */
751          if (hp == NULL)
752                  return(0);
# Line 518 | Line 760 | doambient(                             /* compute ambient component */
760                  dg[0] = dg[1] = 0.0;
761                                          /* sample the hemisphere */
762          acol[0] = acol[1] = acol[2] = 0.0;
763 +        cnt = 0;
764          for (i = hp->ns; i--; )
765                  for (j = hp->ns; j--; )
766                          if ((ap = ambsample(hp, i, j)) != NULL) {
# Line 529 | Line 772 | doambient(                             /* compute ambient component */
772                  free(hp);
773                  return(0);              /* no valid samples */
774          }
775 +        if (cnt < hp->ns*hp->ns) {      /* incomplete sampling? */
776 +                copycolor(rcol, acol);
777 +                free(hp);
778 +                return(-1);             /* return value w/o Hessian */
779 +        }
780 +        cnt = ambssamp*wt + 0.5;        /* perform super-sampling? */
781 +        if (cnt > 0)
782 +                ambsupersamp(acol, hp, cnt);
783          copycolor(rcol, acol);          /* final indirect irradiance/PI */
784 <        if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
534 <                        (ra == NULL) & (pg == NULL) & (dg == NULL)) {
784 >        if ((ra == NULL) & (pg == NULL) & (dg == NULL)) {
785                  free(hp);
786                  return(-1);             /* no radius or gradient calc. */
787          }
788 <        if (bright(acol) > FTINY)       /* normalize Y values */
789 <                d = cnt/bright(acol);
790 <        else
791 <                d = 0.0;
788 >        if ((d = bright(acol)) > FTINY) {       /* normalize Y values */
789 >                d = 0.99*(hp->ns*hp->ns)/d;
790 >                K = 0.01;
791 >        } else {                        /* or fall back on geometric Hessian */
792 >                K = 1.0;
793 >                pg = NULL;
794 >                dg = NULL;
795 >        }
796          ap = hp->sa;                    /* relative Y channel from here on... */
797          for (i = hp->ns*hp->ns; i--; ap++)
798 <                colval(ap->v,CIEY) = bright(ap->v)*d + 0.01;
798 >                colval(ap->v,CIEY) = bright(ap->v)*d + K;
799  
800          if (uv == NULL)                 /* make sure we have axis pointers */
801                  uv = my_uv;
# Line 552 | Line 806 | doambient(                             /* compute ambient component */
806                  ambdirgrad(hp, uv, dg);
807  
808          if (ra != NULL) {               /* scale/clamp radii */
809 +                if (pg != NULL) {
810 +                        if (ra[0]*(d = fabs(pg[0])) > 1.0)
811 +                                ra[0] = 1.0/d;
812 +                        if (ra[1]*(d = fabs(pg[1])) > 1.0)
813 +                                ra[1] = 1.0/d;
814 +                        if (ra[0] > ra[1])
815 +                                ra[0] = ra[1];
816 +                }
817                  if (ra[0] < minarad) {
818                          ra[0] = minarad;
819                          if (ra[1] < minarad)
820                                  ra[1] = minarad;
559 –                                        /* cap gradient if necessary */
560 –                        if (pg != NULL) {
561 –                                d = pg[0]*pg[0]*ra[0]*ra[0] +
562 –                                                pg[1]*pg[1]*ra[1]*ra[1];
563 –                                if (d > 1.0) {
564 –                                        d = 1.0/sqrt(d);
565 –                                        pg[0] *= d;
566 –                                        pg[1] *= d;
567 –                                }
568 –                        }
821                  }
822                  ra[0] *= d = 1.0/sqrt(sqrt(wt));
823                  if ((ra[1] *= d) > 2.0*ra[0])
# Line 574 | Line 826 | doambient(                             /* compute ambient component */
826                          ra[1] = maxarad;
827                          if (ra[0] > maxarad)
828                                  ra[0] = maxarad;
829 +                }
830 +                if (pg != NULL) {       /* cap gradient if necessary */
831 +                        d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1];
832 +                        if (d > 1.0) {
833 +                                d = 1.0/sqrt(d);
834 +                                pg[0] *= d;
835 +                                pg[1] *= d;
836 +                        }
837                  }
838          }
839          free(hp);                       /* clean up and return */

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines