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.45 by greg, Thu May 1 22:34:25 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 */
# Line 34 | Line 61 | typedef struct {
61          AMBSAMP sa[1];          /* sample array (extends struct) */
62   }  AMBHEMI;             /* ambient sample hemisphere */
63  
64 < #define ambsam(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, 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 106 | Line 173 | getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n)
173                  scalecolor(arp->rcoef, 1./AVGREFL);
174          }
175          hlist[0] = hp->rp->rno;
176 <        hlist[1] = i;
177 <        hlist[2] = j;
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))
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))
182 >                if ((spt[1] < 0.1) | (spt[1] >= 0.9))
183                          spt[1] = 0.1 + 0.8*frandom();
184          }
185 <        SDsquare2disk(spt, (i+spt[0])/hp->ns, (j+spt[1])/hp->ns);
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                  arp->rdir[ii] = spt[0]*hp->ux[ii] +
189                                  spt[1]*hp->uy[ii] +
190                                  zd*hp->rp->ron[ii];
191          checknorm(arp->rdir);
192 <        dimlist[ndims++] = i*hp->ns + j + 90171;
192 >        dimlist[ndims++] = ambndx(hp,i,j) + 90171;
193          rayvalue(arp);                  /* evaluate ray */
194          ndims--;                        /* apply coefficient */
195          multcolor(arp->rcol, arp->rcoef);
# Line 245 | Line 312 | ambsupersamp(double acol[3], AMBHEMI *hp, int cnt)
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 <        double  rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2;
397 <        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);
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);
# Line 266 | Line 420 | comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop
420          ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r +
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 293 | Line 448 | comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm)
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);
# Line 316 | 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 338 | 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 356 | Line 516 | comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm)
516          double  f1;
517          int     i;
518  
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*ncp[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 375 | 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 384 | Line 548 | add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F
548   }
549  
550  
387 – /* Return brightness of furthest ambient sample */
388 – static COLORV
389 – back_ambval(AMBSAMP *ap1, AMBSAMP *ap2, AMBSAMP *ap3, FVECT orig)
390 – {
391 –        COLORV  vback;
392 –        FVECT   vec;
393 –        double  d2, d2best;
394 –
395 –        VSUB(vec, ap1->p, orig);
396 –        d2best = DOT(vec,vec);
397 –        vback = colval(ap1->v,CIEY);
398 –        VSUB(vec, ap2->p, orig);
399 –        d2 = DOT(vec,vec);
400 –        if (d2 > d2best) {
401 –                d2best = d2;
402 –                vback = colval(ap2->v,CIEY);
403 –        }
404 –        VSUB(vec, ap3->p, orig);
405 –        d2 = DOT(vec,vec);
406 –        if (d2 > d2best)
407 –                return(colval(ap3->v,CIEY));
408 –        return(vback);
409 – }
410 –
411 –
551   /* Compute anisotropic radii and eigenvector directions */
552   static int
553   eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3])
# Line 469 | 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 490 | 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, ambsam(hp,0,j).p,
496 <                                ambsam(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 503 | 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, ambsam(hp,i,0).p,
507 <                        ambsam(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 512 | 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(&ambsam(hp,i,j), &ambsam(hp,i,j+1),
517 <                                        &ambsam(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, ambsam(hp,i,j+1).p,
520 <                                ambsam(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);
# Line 529 | Line 667 | ambHessian(                            /* anisotropic radii & pos. gradient */
667                      add2gradient(gradient, gradrow[j], graddia, gradcol, backg);
668                  }
669                                          /* initialize edge in next row */
670 <                comp_fftri(&fftr, ambsam(hp,i+1,j+1).p,
533 <                                ambsam(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(&ambsam(hp,i,j+1), &ambsam(hp,i+1,j+1),
677 <                                        &ambsam(hp,i+1,j), hp->rp->rop);
541 <                comp_fftri(&fftr, ambsam(hp,i,j+1).p, ambsam(hp,i+1,j+1).p,
542 <                                hp->rp->rop);
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 559 | 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);

Diff Legend

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