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.21 by greg, Mon Dec 31 18:19:42 2007 UTC vs.
Revision 2.26 by greg, Wed Apr 16 20:32:00 2014 UTC

# Line 10 | Line 10 | static const char      RCSid[] = "$Id$";
10   #include "copyright.h"
11  
12   #include  "ray.h"
13 –
13   #include  "ambient.h"
15 –
14   #include  "random.h"
15  
16 + #ifdef NEWAMB
17  
18 + extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
19 +
20 + typedef struct {
21 +        RAY     *rp;            /* originating ray sample */
22 +        FVECT   ux, uy;         /* tangent axis directions */
23 +        int     ns;             /* number of samples per axis */
24 +        COLOR   acoef;          /* division contribution coefficient */
25 +        struct s_ambsamp {
26 +                COLOR   v;              /* hemisphere sample value */
27 +                float   p[3];           /* intersection point */
28 +        } sa[1];                /* sample array (extends struct) */
29 + }  AMBHEMI;             /* ambient sample hemisphere */
30 +
31 + #define ambsamp(h,i,j)  (h)->sa[(i)*(h)->ns + (j)]
32 +
33 +
34 + static AMBHEMI *
35 + inithemi(                       /* initialize sampling hemisphere */
36 +        COLOR   ac,
37 +        RAY     *r,
38 +        double  wt
39 + )
40 + {
41 +        AMBHEMI *hp;
42 +        double  d;
43 +        int     n, i;
44 +                                        /* set number of divisions */
45 +        if (ambacc <= FTINY &&
46 +                        wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
47 +                wt = d;                 /* avoid ray termination */
48 +        n = sqrt(ambdiv * wt) + 0.5;
49 +        i = 1 + 4*(ambacc > FTINY);     /* minimum number of samples */
50 +        if (n < i)
51 +                n = i;
52 +                                        /* allocate sampling array */
53 +        hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) +
54 +                                sizeof(struct s_ambsamp)*(n*n - 1));
55 +        if (hp == NULL)
56 +                return(NULL);
57 +        hp->rp = r;
58 +        hp->ns = n;
59 +                                        /* assign coefficient */
60 +        copycolor(hp->acoef, ac);
61 +        d = 1.0/(n*n);
62 +        scalecolor(hp->acoef, d);
63 +                                        /* make tangent axes */
64 +        hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
65 +        for (i = 0; i < 3; i++)
66 +                if (r->rn[i] < 0.6 && r->rn[i] > -0.6)
67 +                        break;
68 +        if (i >= 3)
69 +                error(CONSISTENCY, "bad ray direction in inithemi()");
70 +        hp->uy[i] = 1.0;
71 +        VCROSS(hp->ux, hp->uy, r->rn);
72 +        normalize(hp->ux);
73 +        VCROSS(hp->uy, r->rn, hp->ux);
74 +                                        /* we're ready to sample */
75 +        return(hp);
76 + }
77 +
78 +
79 + static int
80 + ambsample(                              /* sample an ambient direction */
81 +        AMBHEMI *hp,
82 +        int     i,
83 +        int     j,
84 + )
85 + {
86 +        struct s_ambsamp        *ap = &ambsamp(hp,i,j);
87 +        RAY                     ar;
88 +        int                     hlist[3];
89 +        double                  spt[2], dz;
90 +        int                     ii;
91 +                                        /* ambient coefficient for weight */
92 +        if (ambacc > FTINY)
93 +                setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
94 +        else
95 +                copycolor(ar.rcoef, hp->acoef);
96 +        if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) {
97 +                setcolor(ap->v, 0., 0., 0.);
98 +                ap->r = 0.;
99 +                return(0);              /* no sample taken */
100 +        }
101 +        if (ambacc > FTINY) {
102 +                multcolor(ar.rcoef, hp->acoef);
103 +                scalecolor(ar.rcoef, 1./AVGREFL);
104 +        }
105 +                                        /* generate hemispherical sample */
106 +        SDsquare2disk(spt, (i+frandom())/hp->ns, (j+frandom())/hp->ns);
107 +        zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
108 +        for (ii = 3; ii--; )
109 +                ar.rdir[ii] =   spt[0]*hp->ux[ii] +
110 +                                spt[1]*hp->uy[ii] +
111 +                                zd*hp->rp->ron[ii];
112 +        checknorm(ar.rdir);
113 +        dimlist[ndims++] = i*hp->ns + j + 90171;
114 +        rayvalue(&ar);                  /* evaluate ray */
115 +        ndims--;
116 +        multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
117 +        copycolor(ap->v, ar.rcol);
118 +        if (ar.rt > 20.0*maxarad)       /* limit vertex distance */
119 +                ar.rt = 20.0*maxarad;
120 +        VSUM(ap->p, ar.rorg, ar.rdir, ar.rt);
121 +        return(1);
122 + }
123 +
124 +
125 + static void
126 + ambHessian(                             /* anisotropic radii & pos. gradient */
127 +        AMBHEMI *hp,
128 +        FVECT   uv[2],                  /* returned */
129 +        float   ra[2],                  /* returned */
130 +        float   pg[2]                   /* returned */
131 + )
132 + {
133 +        if (ra != NULL) {               /* compute Hessian-derived radii */
134 +        } else {                        /* else copy original tangent axes */
135 +                VCOPY(uv[0], hp->ux);
136 +                VCOPY(uv[1], hp->uy);
137 +        }
138 +        if (pg == NULL)                 /* no position gradient requested? */
139 +                return;
140 + }
141 +
142 + int
143 + doambient(                              /* compute ambient component */
144 +        COLOR   rcol,                   /* input/output color */
145 +        RAY     *r,
146 +        double  wt,
147 +        FVECT   uv[2],                  /* returned */
148 +        float   ra[2],                  /* returned */
149 +        float   pg[2],                  /* returned */
150 +        float   dg[2]                   /* returned */
151 + )
152 + {
153 +        int                     cnt = 0;
154 +        FVECT                   my_uv[2];
155 +        AMBHEMI                 *hp;
156 +        double                  d, acol[3];
157 +        struct s_ambsamp        *ap;
158 +        int                     i, j;
159 +                                        /* initialize */
160 +        if ((hp = inithemi(rcol, r, wt)) == NULL)
161 +                return(0);
162 +        if (uv != NULL)
163 +                memset(uv, 0, sizeof(FVECT)*2);
164 +        if (ra != NULL)
165 +                ra[0] = ra[1] = 0.0;
166 +        if (pg != NULL)
167 +                pg[0] = pg[1] = 0.0;
168 +        if (dg != NULL)
169 +                dg[0] = dg[1] = 0.0;
170 +                                        /* sample the hemisphere */
171 +        acol[0] = acol[1] = acol[2] = 0.0;
172 +        for (i = hemi.ns; i--; )
173 +                for (j = hemi.ns; j--; )
174 +                        if (ambsample(hp, i, j)) {
175 +                                ap = &ambsamp(hp,i,j);
176 +                                addcolor(acol, ap->v);
177 +                                ++cnt;
178 +                        }
179 +        if (!cnt) {
180 +                setcolor(rcol, 0.0, 0.0, 0.0);
181 +                free(hp);
182 +                return(0);              /* no valid samples */
183 +        }
184 +        d = 1.0 / cnt;                  /* final indirect irradiance/PI */
185 +        acol[0] *= d; acol[1] *= d; acol[2] *= d;
186 +        copycolor(rcol, acol);
187 +        if (cnt < hp->ns*hp->ns ||      /* incomplete sampling? */
188 +                        (ra == NULL) & (pg == NULL) & (dg == NULL)) {
189 +                free(hp);
190 +                return(-1);             /* no radius or gradient calc. */
191 +        }
192 +        d = 0.01 * bright(rcol);        /* add in 1% before Hessian comp. */
193 +        if (d < FTINY) d = FTINY;
194 +        ap = hp->sa;                    /* using Y channel from here on... */
195 +        for (i = hp->ns*hp->ns; i--; ap++)
196 +                colval(ap->v,CIEY) = bright(ap->v) + d;
197 +
198 +        if (uv == NULL)                 /* make sure we have axis pointers */
199 +                uv = my_uv;
200 +                                        /* compute radii & pos. gradient */
201 +        ambHessian(hp, uv, ra, pg);
202 +        if (dg != NULL)                 /* compute direction gradient */
203 +                ambdirgrad(hp, uv, dg);
204 +        if (ra != NULL) {               /* adjust/clamp radii */
205 +                d = pow(wt, -0.25);
206 +                if ((ra[0] *= d) > maxarad)
207 +                        ra[0] = maxarad;
208 +                if ((ra[1] *= d) > 2.0*ra[0])
209 +                        ra[1] = 2.0*ra[0];
210 +        }
211 +        free(hp);                       /* clean up and return */
212 +        return(1);
213 + }
214 +
215 +
216 + #else /* ! NEWAMB */
217 +
218 +
219   void
220   inithemi(                       /* initialize sampling hemisphere */
221 <        register AMBHEMI  *hp,
221 >        AMBHEMI  *hp,
222          COLOR ac,
223          RAY  *r,
224          double  wt
225   )
226   {
227          double  d;
228 <        register int  i;
228 >        int  i;
229                                          /* set number of divisions */
230          if (ambacc <= FTINY &&
231                          wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
# Line 58 | Line 258 | inithemi(                      /* initialize sampling hemisphere */
258  
259   int
260   divsample(                              /* sample a division */
261 <        register AMBSAMP  *dp,
261 >        AMBSAMP  *dp,
262          AMBHEMI  *h,
263          RAY  *r
264   )
# Line 69 | Line 269 | divsample(                             /* sample a division */
269          double  xd, yd, zd;
270          double  b2;
271          double  phi;
272 <        register int  i;
272 >        int  i;
273                                          /* ambient coefficient for weight */
274          if (ambacc > FTINY)
275                  setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
# Line 94 | Line 294 | divsample(                             /* sample a division */
294                  ar.rdir[i] =    xd*h->ux[i] +
295                                  yd*h->uy[i] +
296                                  zd*h->uz[i];
297 +        checknorm(ar.rdir);
298          dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
299          rayvalue(&ar);
300          ndims--;
# Line 138 | Line 339 | ambnorm(                               /* standard order */
339   {
340          const AMBSAMP   *d1 = (const AMBSAMP *)p1;
341          const AMBSAMP   *d2 = (const AMBSAMP *)p2;
342 <        register int    c;
342 >        int     c;
343  
344          if ( (c = d1->t - d2->t) )
345                  return(c);
# Line 148 | Line 349 | ambnorm(                               /* standard order */
349  
350   double
351   doambient(                              /* compute ambient component */
352 <        COLOR  acol,
352 >        COLOR  rcol,
353          RAY  *r,
354          double  wt,
355          FVECT  pg,
356          FVECT  dg
357   )
358   {
359 <        double  b, d;
359 >        double  b, d=0;
360          AMBHEMI  hemi;
361          AMBSAMP  *div;
362          AMBSAMP  dnew;
363 <        register AMBSAMP  *dp;
363 >        double  acol[3];
364 >        AMBSAMP  *dp;
365          double  arad;
366          int  divcnt;
367 <        register int  i, j;
367 >        int  i, j;
368                                          /* initialize hemisphere */
369 <        inithemi(&hemi, acol, r, wt);
369 >        inithemi(&hemi, rcol, r, wt);
370          divcnt = hemi.nt * hemi.np;
371                                          /* initialize */
372          if (pg != NULL)
373                  pg[0] = pg[1] = pg[2] = 0.0;
374          if (dg != NULL)
375                  dg[0] = dg[1] = dg[2] = 0.0;
376 <        setcolor(acol, 0.0, 0.0, 0.0);
376 >        setcolor(rcol, 0.0, 0.0, 0.0);
377          if (divcnt == 0)
378                  return(0.0);
379                                          /* allocate super-samples */
# Line 183 | Line 385 | doambient(                             /* compute ambient component */
385                  div = NULL;
386                                          /* sample the divisions */
387          arad = 0.0;
388 +        acol[0] = acol[1] = acol[2] = 0.0;
389          if ((dp = div) == NULL)
390                  dp = &dnew;
391          divcnt = 0;
# Line 264 | Line 467 | doambient(                             /* compute ambient component */
467                  }
468                  free((void *)div);
469          }
470 +        copycolor(rcol, acol);
471          if (arad <= FTINY)
472                  arad = maxarad;
473          else
# Line 290 | Line 494 | doambient(                             /* compute ambient component */
494   void
495   comperrs(                       /* compute initial error estimates */
496          AMBSAMP  *da,   /* assumes standard ordering */
497 <        register AMBHEMI  *hp
497 >        AMBHEMI  *hp
498   )
499   {
500          double  b, b2;
501          int  i, j;
502 <        register AMBSAMP  *dp;
502 >        AMBSAMP  *dp;
503                                  /* sum differences from neighbors */
504          dp = da;
505          for (i = 0; i < hp->nt; i++)
# Line 343 | Line 547 | void
547   posgradient(                                    /* compute position gradient */
548          FVECT  gv,
549          AMBSAMP  *da,                   /* assumes standard ordering */
550 <        register AMBHEMI  *hp
550 >        AMBHEMI  *hp
551   )
552   {
553 <        register int  i, j;
553 >        int  i, j;
554          double  nextsine, lastsine, b, d;
555          double  mag0, mag1;
556          double  phi, cosp, sinp, xd, yd;
557 <        register AMBSAMP  *dp;
557 >        AMBSAMP  *dp;
558  
559          xd = yd = 0.0;
560          for (j = 0; j < hp->np; j++) {
# Line 401 | Line 605 | void
605   dirgradient(                                    /* compute direction gradient */
606          FVECT  gv,
607          AMBSAMP  *da,                   /* assumes standard ordering */
608 <        register AMBHEMI  *hp
608 >        AMBHEMI  *hp
609   )
610   {
611 <        register int  i, j;
611 >        int  i, j;
612          double  mag;
613          double  phi, xd, yd;
614 <        register AMBSAMP  *dp;
614 >        AMBSAMP  *dp;
615  
616          xd = yd = 0.0;
617          for (j = 0; j < hp->np; j++) {
# Line 430 | Line 634 | dirgradient(                                   /* compute direction gradient */
634          for (i = 0; i < 3; i++)
635                  gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
636   }
637 +
638 + #endif  /* ! NEWAMB */

Diff Legend

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