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

Comparing ray/src/gen/mksource.c (file contents):
Revision 2.4 by greg, Fri Sep 23 19:22:37 2005 UTC vs.
Revision 2.5 by greg, Fri Aug 24 04:41:36 2007 UTC

# Line 17 | Line 17 | static const char RCSid[] = "$Id$";
17   /* Data structure for geodesic samples */
18  
19   typedef struct tritree {
20 <        FVECT           gdv[3];         /* spherical triangle vertex direc. */
21 <        FVECT           sd;             /* sample direction if leaf */
20 >        int32           gdv[3];         /* spherical triangle vertex direc. */
21 >        int32           sd;             /* sample direction if leaf */
22          struct tritree  *kid;           /* 4 children if branch node */
23          COLR            val;            /* sampled color value */
24   } TRITREE;
25  
26   typedef struct lostlight {
27          struct lostlight        *next;  /* next in list */
28 <        FVECT           sd;             /* lost source direction */
28 >        int32           sd;             /* lost source direction */
29          COLOR           intens;         /* output times solid angle */
30   } LOSTLIGHT;
31  
# Line 56 | Line 56 | vol_sign(const FVECT v1, const FVECT v2, const FVECT v
56  
57   /* Is the given direction contained within the specified spherical triangle? */
58   int
59 < intriv(FVECT tri[3], const FVECT sdir)
59 > intriv(const int32 trid[3], const FVECT sdir)
60   {
61          int     sv[3];
62 +        FVECT   tri[3];
63  
64 +        decodedir(tri[0], trid[0]);
65 +        decodedir(tri[1], trid[1]);
66 +        decodedir(tri[2], trid[2]);
67          sv[0] = vol_sign(sdir, tri[0], tri[1]);
68          sv[1] = vol_sign(sdir, tri[1], tri[2]);
69          sv[2] = vol_sign(sdir, tri[2], tri[0]);
# Line 91 | Line 95 | leafsample(TRITREE *leaf)
95   {
96          RAY     myray;
97          RREAL   wt[3];
98 +        FVECT   sdir;
99          int     i, j;
100                                                  /* random point on triangle */
101          i = random() % 3;
# Line 98 | Line 103 | leafsample(TRITREE *leaf)
103          j = random() & 1;
104          wt[(i+2-j)%3] = 1. - wt[i] -
105                          (wt[(i+1+j)%3] = (1.-wt[i])*frandom());
106 <        leaf->sd[0] = leaf->sd[1] = leaf->sd[2] = .0;
107 <        for (i = 0; i < 3; i++)
108 <                VSUM(leaf->sd, leaf->sd, leaf->gdv[i], wt[i]);
109 <        normalize(leaf->sd);                    /* record sample direction */
106 >        sdir[0] = sdir[1] = sdir[2] = .0;
107 >        for (i = 0; i < 3; i++) {
108 >                FVECT   vt;
109 >                decodedir(vt, leaf->gdv[i]);
110 >                VSUM(sdir, sdir, vt, wt[i]);
111 >        }
112 >        normalize(sdir);                        /* record sample direction */
113 >        leaf->sd = encodedir(sdir);
114                                                  /* evaluate at inf. */
115 <        VSUM(myray.rorg, scene_cent, leaf->sd, scene_rad);
116 <        VCOPY(myray.rdir, leaf->sd);
115 >        VSUM(myray.rorg, scene_cent, sdir, scene_rad);
116 >        VCOPY(myray.rdir, sdir);
117          myray.rmax = 0.;
118          ray_trace(&myray);
119          setcolr(leaf->val, colval(myray.rcol,RED),
# Line 114 | Line 123 | leafsample(TRITREE *leaf)
123  
124   /* Initialize a branch node contained in the given spherical triangle */
125   void
126 < subdivide(TRITREE *branch, FVECT dv[3])
126 > subdivide(TRITREE *branch, const int32 dv[3])
127   {
128 <        FVECT   sdv[3];
128 >        FVECT   dvv[3], sdv[3];
129 >        int32   sd[3];
130          int     i;
131          
132 <        for (i = 0; i < 3; i++)         /* copy spherical triangle */
133 <                VCOPY(branch->gdv[i], dv[i]);
132 >        for (i = 0; i < 3; i++) {       /* copy spherical triangle */
133 >                branch->gdv[i] = dv[i];
134 >                decodedir(dvv[i], dv[i]);
135 >        }
136          for (i = 0; i < 3; i++) {       /* create new vertices */
137                  int     j = (i+1)%3;
138 <                VADD(sdv[i], dv[i], dv[j]);
138 >                VADD(sdv[i], dvv[i], dvv[j]);
139                  normalize(sdv[i]);
140 +                sd[i] = encodedir(sdv[i]);
141          }
142                                          /* allocate leaves */
143          branch->kid = (TRITREE *)calloc(4, sizeof(TRITREE));
144          if (branch->kid == NULL)
145                  error(SYSTEM, "out of memory in subdivide()");
146                                          /* assign subtriangle directions */
147 <        VCOPY(branch->kid[0].gdv[0], dv[0]);
148 <        VCOPY(branch->kid[0].gdv[1], sdv[0]);
149 <        VCOPY(branch->kid[0].gdv[2], sdv[2]);
150 <        VCOPY(branch->kid[1].gdv[0], sdv[0]);
151 <        VCOPY(branch->kid[1].gdv[1], dv[1]);
152 <        VCOPY(branch->kid[1].gdv[2], sdv[1]);
153 <        VCOPY(branch->kid[2].gdv[0], sdv[1]);
154 <        VCOPY(branch->kid[2].gdv[1], dv[2]);
155 <        VCOPY(branch->kid[2].gdv[2], sdv[2]);
156 <        VCOPY(branch->kid[3].gdv[0], sdv[0]);
157 <        VCOPY(branch->kid[3].gdv[1], sdv[1]);
158 <        VCOPY(branch->kid[3].gdv[2], sdv[2]);
147 >        branch->kid[0].gdv[0] = dv[0];
148 >        branch->kid[0].gdv[1] = sd[0];
149 >        branch->kid[0].gdv[2] = sd[2];
150 >        branch->kid[1].gdv[0] = sd[0];
151 >        branch->kid[1].gdv[1] = dv[1];
152 >        branch->kid[1].gdv[2] = sd[1];
153 >        branch->kid[2].gdv[0] = sd[1];
154 >        branch->kid[2].gdv[1] = dv[2];
155 >        branch->kid[2].gdv[2] = sd[2];
156 >        branch->kid[3].gdv[0] = sd[0];
157 >        branch->kid[3].gdv[1] = sd[1];
158 >        branch->kid[3].gdv[2] = sd[2];
159   }
160  
161   /* Recursively subdivide the given node to the specified quadtree depth */
# Line 155 | Line 168 | branchsample(TRITREE *node, int depth)
168                  return;
169          if (isleaf(node)) {                     /* subdivide leaf node */
170                  TRITREE branch, *moved_leaf;
171 +                FVECT   sdir;
172                  subdivide(&branch, node->gdv);
173 <                moved_leaf = findleaf(&branch, node->sd);
173 >                decodedir(sdir, node->sd);
174 >                moved_leaf = findleaf(&branch, sdir);
175                  if (moved_leaf != NULL) {       /* bequeath old sample */
176 <                        VCOPY(moved_leaf->sd, node->sd);
176 >                        moved_leaf->sd = node->sd;
177                          copycolr(moved_leaf->val, node->val);
178                  }
179                  for (i = 0; i < 4; i++)         /* compute new samples */
# Line 193 | Line 208 | geosample(int nsamps)
208          trunk[1][1] = sqrt(1. - trunk[1][2]*trunk[1][2]);
209          spinvector(trunk[2], trunk[1], trunk[0], 2.*PI/3.);
210          spinvector(trunk[3], trunk[1], trunk[0], 4.*PI/3.);
211 <        VCOPY(tree->gdv[0], trunk[0]);
197 <        VCOPY(tree->gdv[1], trunk[0]);
198 <        VCOPY(tree->gdv[2], trunk[0]);
211 >        tree->gdv[0] = tree->gdv[1] = tree->gdv[2] = encodedir(trunk[0]);
212          tree->kid = (TRITREE *)calloc(NTRUNKBR, sizeof(TRITREE));
213          if (tree->kid == NULL) goto memerr;
214                                          /* grow our tree from trunk */
215          for (i = 0; i < NTRUNKBR; i++) {
216                  for (j = 0; j < 3; j++)         /* XXX works for tetra only */
217 <                        VCOPY(tree->kid[i].gdv[j], trunk[(i+j)%NTRUNKVERT]);
217 >                        tree->kid[i].gdv[j] = encodedir(trunk[(i+j)%NTRUNKVERT]);
218                  leafsample(&tree->kid[i]);
219                  branchsample(&tree->kid[i], depth);
220          }
# Line 269 | Line 282 | findemax(TRITREE *node, int *expp)
282  
283   /* Compute solid angle of spherical triangle (approx.) */
284   double
285 < tri_omegav(FVECT v[3])
285 > tri_omegav(const int32 vd[3])
286   {
287 <        FVECT   e1, e2, vcross;
288 <        
287 >        FVECT   v[3], e1, e2, vcross;
288 >
289 >        decodedir(v[0], vd[0]);
290 >        decodedir(v[1], vd[1]);
291 >        decodedir(v[2], vd[2]);
292          VSUB(e1, v[1], v[0]);
293          VSUB(e2, v[2], v[1]);
294          fcross(vcross, e1, e2);
295          return(.5*VLEN(vcross));
296   }
297  
298 < /* Sum intensity times direction for non-zero leaves */
298 > /* Sum intensity times direction for above-threshold perimiter within radius */
299   void
300   vector_sum(FVECT vsum, TRITREE *node,
301 <                const FVECT cent, double mincos, int ethresh)
301 >                FVECT cent, double maxr2, int ethresh)
302   {
303          if (isleaf(node)) {
304                  double  intens;
305 +                FVECT   sdir;
306                  if (node->val[EXP] < ethresh)
307 <                        return;
308 <                if (DOT(node->sd,cent) < mincos)
309 <                        return;
307 >                        return;                         /* below threshold */
308 >                if (fdir2diff(node->sd,cent) > maxr2)
309 >                        return;                         /* too far away */
310                  intens = colrval(node->val,GRN) * tri_omegav(node->gdv);
311 <                VSUM(vsum, vsum, node->sd, intens);
311 >                decodedir(sdir, node->sd);
312 >                VSUM(vsum, vsum, sdir, intens);
313                  return;
314          }
315 <        if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
316 <                        DOT(node->gdv[0],cent) > mincos &&
317 <                        DOT(node->gdv[1],cent) > mincos &&
318 <                        DOT(node->gdv[2],cent) > mincos)
319 <                return;
320 <        vector_sum(vsum, &node->kid[0], cent, mincos, ethresh);
321 <        vector_sum(vsum, &node->kid[1], cent, mincos, ethresh);
322 <        vector_sum(vsum, &node->kid[2], cent, mincos, ethresh);
323 <        vector_sum(vsum, &node->kid[3], cent, mincos, ethresh);
315 >        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
316 >                        fdir2diff(node->gdv[0],cent) < maxr2 &&
317 >                        fdir2diff(node->gdv[1],cent) < maxr2 &&
318 >                        fdir2diff(node->gdv[2],cent) < maxr2)
319 >                return;                                 /* containing node */
320 >        vector_sum(vsum, &node->kid[0], cent, maxr2, ethresh);
321 >        vector_sum(vsum, &node->kid[1], cent, maxr2, ethresh);
322 >        vector_sum(vsum, &node->kid[2], cent, maxr2, ethresh);
323 >        vector_sum(vsum, &node->kid[3], cent, maxr2, ethresh);
324   }
325  
326   /* Claim source contributions within the given solid angle */
327   void
328 < claimlight(COLOR intens, TRITREE *node, const FVECT cent, double mincos)
328 > claimlight(COLOR intens, TRITREE *node, FVECT cent, double maxr2)
329   {
330          int     remaining;
331          int     i;
332          if (isleaf(node)) {             /* claim contribution */
333                  COLOR   contrib;
334                  if (node->val[EXP] <= 0)
335 <                        return;
336 <                if (DOT(node->sd,cent) < mincos)
337 <                        return;
335 >                        return;         /* already claimed */
336 >                if (fdir2diff(node->sd,cent) > maxr2)
337 >                        return;         /* too far away */
338                  colr_color(contrib, node->val);
339                  scalecolor(contrib, tri_omegav(node->gdv));
340                  addcolor(intens, contrib);
341                  copycolr(node->val, blkclr);
342                  return;
343          }
344 <        if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
345 <                        DOT(node->gdv[0],cent) > mincos &&
346 <                        DOT(node->gdv[1],cent) > mincos &&
347 <                        DOT(node->gdv[2],cent) > mincos)
348 <                return;
344 >        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
345 >                        fdir2diff(node->gdv[0],cent) < maxr2 &&
346 >                        fdir2diff(node->gdv[1],cent) < maxr2 &&
347 >                        fdir2diff(node->gdv[2],cent) < maxr2)
348 >                return;                 /* previously claimed node */
349          remaining = 0;                  /* recurse on children */
350          for (i = 0; i < 4; i++) {
351 <                claimlight(intens, &node->kid[i], cent, mincos);
351 >                claimlight(intens, &node->kid[i], cent, maxr2);
352                  if (!isleaf(&node->kid[i]) || node->kid[i].val[EXP] != 0)
353                          ++remaining;
354          }
# Line 339 | Line 357 | claimlight(COLOR intens, TRITREE *node, const FVECT ce
357                                          /* consolidate empties */
358          free((void *)node->kid); node->kid = NULL;
359          copycolr(node->val, blkclr);
360 <        VCOPY(node->sd, node->gdv[0]);  /* doesn't really matter */
360 >        node->sd = node->gdv[0];        /* doesn't really matter */
361   }
362  
363   /* Add lost light contribution to the given list */
364   void
365 < add2lost(LOSTLIGHT **llp, COLOR intens, const FVECT cent)
365 > add2lost(LOSTLIGHT **llp, COLOR intens, FVECT cent)
366   {
367          LOSTLIGHT       *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
368  
369          if (newll == NULL)
370                  return;
371          copycolor(newll->intens, intens);
372 <        VCOPY(newll->sd, cent);
372 >        newll->sd = encodedir(cent);
373          newll->next = *llp;
374          *llp = newll;
375   }
376  
377   /* Check lost light list for contributions */
378   void
379 < getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cent, double omega)
379 > getlost(LOSTLIGHT **llp, COLOR intens, FVECT cent, double omega)
380   {
381 <        const double    mincos = 1. - omega/(2.*PI);
381 >        const double    maxr2 = omega/PI;
382          LOSTLIGHT       lhead, *lastp, *thisp;
383  
384          lhead.next = *llp;
385          lastp = &lhead;
386          while ((thisp = lastp->next) != NULL)
387 <                if (DOT(thisp->sd,cent) >= mincos) {
387 >                if (fdir2diff(thisp->sd,cent) <= maxr2) {
388                          LOSTLIGHT       *mynext = thisp->next;
389                          addcolor(intens, thisp->intens);
390                          free((void *)thisp);
# Line 419 | Line 437 | mksources(TRITREE *samptree, double thresh, double max
437                  if (startleaf == NULL)
438                          break;
439                                          /* claim it */
440 <                VCOPY(curcent, startleaf->sd);
440 >                decodedir(curcent, startleaf->sd);
441                  curomega = tri_omegav(startleaf->gdv);
442                  currad = sqrt(curomega/PI);
443                  growstep = 3.*currad;
# Line 435 | Line 453 | mksources(TRITREE *samptree, double thresh, double max
453                          vsum[0] = vsum[1] = vsum[2] = .0;
454                          for (i = 0; i < NTRUNKBR; i++)
455                                  vector_sum(vsum, &samptree->kid[i],
456 <                                                curcent, cos(currad+growstep),
456 >                                        curcent, 2.-2.*cos(currad+growstep),
457                                                  thisethresh);
458                          if (normalize(vsum) == .0)
459                                  break;
# Line 451 | Line 469 | mksources(TRITREE *samptree, double thresh, double max
469                          curomega = 2.*PI*(1. - cos(currad));
470                          for (i = 0; i < NTRUNKBR; i++)
471                                  claimlight(curintens, &samptree->kid[i],
472 <                                                curcent, cos(currad));
472 >                                                curcent, 2.-2.*cos(currad));
473                  } while (curomega < maxomega &&
474                                  bright(curintens)/curomega > thisthresh);
475                  if (bright(curintens) < minintens) {

Diff Legend

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