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.13 by greg, Tue Jun 3 21:31:51 2025 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7  
8   #include "ray.h"
9   #include "random.h"
10 + #include "paths.h"
11   #include "resolu.h"
12  
13   #define NTRUNKBR        4               /* number of branches at trunk */
# Line 17 | Line 18 | static const char RCSid[] = "$Id$";
18   /* Data structure for geodesic samples */
19  
20   typedef struct tritree {
21 <        FVECT           gdv[3];         /* spherical triangle vertex direc. */
22 <        FVECT           sd;             /* sample direction if leaf */
21 >        int32           gdv[3];         /* spherical triangle vertex direc. */
22 >        int32           sd;             /* sample direction if leaf */
23          struct tritree  *kid;           /* 4 children if branch node */
24          COLR            val;            /* sampled color value */
25   } TRITREE;
26  
27   typedef struct lostlight {
28          struct lostlight        *next;  /* next in list */
29 <        FVECT           sd;             /* lost source direction */
29 >        int32           sd;             /* lost source direction */
30          COLOR           intens;         /* output times solid angle */
31   } LOSTLIGHT;
32  
32 extern char     *progname;
33
33   FVECT   scene_cent;             /* center of octree cube */
34   RREAL   scene_rad;              /* radius to get outside cube from center */
35  
# Line 56 | Line 55 | vol_sign(const FVECT v1, const FVECT v2, const FVECT v
55  
56   /* Is the given direction contained within the specified spherical triangle? */
57   int
58 < intriv(FVECT tri[3], const FVECT sdir)
58 > intriv(const int32 trid[3], const FVECT sdir)
59   {
60          int     sv[3];
61 +        FVECT   tri[3];
62  
63 +        decodedir(tri[0], trid[0]);
64 +        decodedir(tri[1], trid[1]);
65 +        decodedir(tri[2], trid[2]);
66          sv[0] = vol_sign(sdir, tri[0], tri[1]);
67          sv[1] = vol_sign(sdir, tri[1], tri[2]);
68          sv[2] = vol_sign(sdir, tri[2], tri[0]);
# Line 91 | Line 94 | leafsample(TRITREE *leaf)
94   {
95          RAY     myray;
96          RREAL   wt[3];
97 +        FVECT   sdir;
98          int     i, j;
99                                                  /* random point on triangle */
100          i = random() % 3;
# Line 98 | Line 102 | leafsample(TRITREE *leaf)
102          j = random() & 1;
103          wt[(i+2-j)%3] = 1. - wt[i] -
104                          (wt[(i+1+j)%3] = (1.-wt[i])*frandom());
105 <        leaf->sd[0] = leaf->sd[1] = leaf->sd[2] = .0;
106 <        for (i = 0; i < 3; i++)
107 <                VSUM(leaf->sd, leaf->sd, leaf->gdv[i], wt[i]);
108 <        normalize(leaf->sd);                    /* record sample direction */
105 >        sdir[0] = sdir[1] = sdir[2] = .0;
106 >        for (i = 0; i < 3; i++) {
107 >                FVECT   vt;
108 >                decodedir(vt, leaf->gdv[i]);
109 >                VSUM(sdir, sdir, vt, wt[i]);
110 >        }
111 >        normalize(sdir);                        /* record sample direction */
112 >        leaf->sd = encodedir(sdir);
113                                                  /* evaluate at inf. */
114 <        VSUM(myray.rorg, scene_cent, leaf->sd, scene_rad);
115 <        VCOPY(myray.rdir, leaf->sd);
114 >        VSUM(myray.rorg, scene_cent, sdir, scene_rad);
115 >        VCOPY(myray.rdir, sdir);
116          myray.rmax = 0.;
117          ray_trace(&myray);
118 <        setcolr(leaf->val, colval(myray.rcol,RED),
111 <                        colval(myray.rcol,GRN),
112 <                        colval(myray.rcol,BLU));
118 >        scolor_colr(leaf->val, myray.rcol);
119   }
120  
121   /* Initialize a branch node contained in the given spherical triangle */
122   void
123 < subdivide(TRITREE *branch, FVECT dv[3])
123 > subdivide(TRITREE *branch, const int32 dv[3])
124   {
125 <        FVECT   sdv[3];
125 >        FVECT   dvv[3], sdv[3];
126 >        int32   sd[3];
127          int     i;
128          
129 <        for (i = 0; i < 3; i++)         /* copy spherical triangle */
130 <                VCOPY(branch->gdv[i], dv[i]);
129 >        for (i = 0; i < 3; i++) {       /* copy spherical triangle */
130 >                branch->gdv[i] = dv[i];
131 >                decodedir(dvv[i], dv[i]);
132 >        }
133          for (i = 0; i < 3; i++) {       /* create new vertices */
134                  int     j = (i+1)%3;
135 <                VADD(sdv[i], dv[i], dv[j]);
135 >                VADD(sdv[i], dvv[i], dvv[j]);
136                  normalize(sdv[i]);
137 +                sd[i] = encodedir(sdv[i]);
138          }
139                                          /* allocate leaves */
140          branch->kid = (TRITREE *)calloc(4, sizeof(TRITREE));
141          if (branch->kid == NULL)
142                  error(SYSTEM, "out of memory in subdivide()");
143                                          /* assign subtriangle directions */
144 <        VCOPY(branch->kid[0].gdv[0], dv[0]);
145 <        VCOPY(branch->kid[0].gdv[1], sdv[0]);
146 <        VCOPY(branch->kid[0].gdv[2], sdv[2]);
147 <        VCOPY(branch->kid[1].gdv[0], sdv[0]);
148 <        VCOPY(branch->kid[1].gdv[1], dv[1]);
149 <        VCOPY(branch->kid[1].gdv[2], sdv[1]);
150 <        VCOPY(branch->kid[2].gdv[0], sdv[1]);
151 <        VCOPY(branch->kid[2].gdv[1], dv[2]);
152 <        VCOPY(branch->kid[2].gdv[2], sdv[2]);
153 <        VCOPY(branch->kid[3].gdv[0], sdv[0]);
154 <        VCOPY(branch->kid[3].gdv[1], sdv[1]);
155 <        VCOPY(branch->kid[3].gdv[2], sdv[2]);
144 >        branch->kid[0].gdv[0] = dv[0];
145 >        branch->kid[0].gdv[1] = sd[0];
146 >        branch->kid[0].gdv[2] = sd[2];
147 >        branch->kid[1].gdv[0] = sd[0];
148 >        branch->kid[1].gdv[1] = dv[1];
149 >        branch->kid[1].gdv[2] = sd[1];
150 >        branch->kid[2].gdv[0] = sd[1];
151 >        branch->kid[2].gdv[1] = dv[2];
152 >        branch->kid[2].gdv[2] = sd[2];
153 >        branch->kid[3].gdv[0] = sd[0];
154 >        branch->kid[3].gdv[1] = sd[1];
155 >        branch->kid[3].gdv[2] = sd[2];
156   }
157  
158   /* Recursively subdivide the given node to the specified quadtree depth */
# Line 155 | Line 165 | branchsample(TRITREE *node, int depth)
165                  return;
166          if (isleaf(node)) {                     /* subdivide leaf node */
167                  TRITREE branch, *moved_leaf;
168 +                FVECT   sdir;
169                  subdivide(&branch, node->gdv);
170 <                moved_leaf = findleaf(&branch, node->sd);
170 >                decodedir(sdir, node->sd);
171 >                moved_leaf = findleaf(&branch, sdir);
172                  if (moved_leaf != NULL) {       /* bequeath old sample */
173 <                        VCOPY(moved_leaf->sd, node->sd);
173 >                        moved_leaf->sd = node->sd;
174                          copycolr(moved_leaf->val, node->val);
175                  }
176                  for (i = 0; i < 4; i++)         /* compute new samples */
# Line 172 | Line 184 | branchsample(TRITREE *node, int depth)
184  
185   /* Sample sphere using triangular geodesic mesh */
186   TRITREE *
187 < geosample(int nsamps)
187 > geosample(long nsamps)
188   {
189          int     depth;
190          TRITREE *tree;
# Line 181 | Line 193 | geosample(int nsamps)
193                                          /* figure out depth */
194          if ((nsamps -= 4) < 0)
195                  error(USER, "minimum number of samples is 4");
196 <        nsamps = nsamps*3/NTRUNKBR;     /* round up */
196 >        nsamps = nsamps*(NTRUNKBR-1)/NTRUNKBR;  /* round up */
197          for (depth = 0; nsamps > 1; depth++)
198                  nsamps >>= 2;
199                                          /* make base tetrahedron */
# Line 193 | Line 205 | geosample(int nsamps)
205          trunk[1][1] = sqrt(1. - trunk[1][2]*trunk[1][2]);
206          spinvector(trunk[2], trunk[1], trunk[0], 2.*PI/3.);
207          spinvector(trunk[3], trunk[1], trunk[0], 4.*PI/3.);
208 <        VCOPY(tree->gdv[0], trunk[0]);
197 <        VCOPY(tree->gdv[1], trunk[0]);
198 <        VCOPY(tree->gdv[2], trunk[0]);
208 >        tree->gdv[0] = tree->gdv[1] = tree->gdv[2] = encodedir(trunk[0]);
209          tree->kid = (TRITREE *)calloc(NTRUNKBR, sizeof(TRITREE));
210          if (tree->kid == NULL) goto memerr;
211                                          /* grow our tree from trunk */
212          for (i = 0; i < NTRUNKBR; i++) {
213                  for (j = 0; j < 3; j++)         /* XXX works for tetra only */
214 <                        VCOPY(tree->kid[i].gdv[j], trunk[(i+j)%NTRUNKVERT]);
214 >                        tree->kid[i].gdv[j] = encodedir(trunk[(i+j)%NTRUNKVERT]);
215                  leafsample(&tree->kid[i]);
216                  branchsample(&tree->kid[i], depth);
217          }
# Line 229 | Line 239 | get_ehisto(const TRITREE *node, long exphisto[256])
239   double
240   get_threshold(const TRITREE *tree)
241   {
242 +        long    samptotal = 0;
243          long    exphisto[256];
233        long    samptotal;
244          int     i;
245                                                  /* compute sample histogram */
246 <        memset((void *)exphisto, 0, sizeof(exphisto));
246 >        memset(exphisto, 0, sizeof(exphisto));
247          for (i = 0; i < NTRUNKBR; i++)
248                  get_ehisto(&tree->kid[i], exphisto);
249                                                  /* use 98th percentile */
# Line 269 | Line 279 | findemax(TRITREE *node, int *expp)
279  
280   /* Compute solid angle of spherical triangle (approx.) */
281   double
282 < tri_omegav(FVECT v[3])
282 > tri_omegav(const int32 vd[3])
283   {
284 <        FVECT   e1, e2, vcross;
285 <        
284 >        FVECT   v[3], e1, e2, vcross;
285 >
286 >        decodedir(v[0], vd[0]);
287 >        decodedir(v[1], vd[1]);
288 >        decodedir(v[2], vd[2]);
289          VSUB(e1, v[1], v[0]);
290          VSUB(e2, v[2], v[1]);
291          fcross(vcross, e1, e2);
292          return(.5*VLEN(vcross));
293   }
294  
295 < /* Sum intensity times direction for non-zero leaves */
295 > /* Sum intensity times direction for above-threshold perimiter within radius */
296   void
297   vector_sum(FVECT vsum, TRITREE *node,
298 <                const FVECT cent, double mincos, int ethresh)
298 >                FVECT cent, double maxr2, int ethresh)
299   {
300          if (isleaf(node)) {
301                  double  intens;
302 +                FVECT   sdir;
303                  if (node->val[EXP] < ethresh)
304 <                        return;
305 <                if (DOT(node->sd,cent) < mincos)
306 <                        return;
304 >                        return;                         /* below threshold */
305 >                if (fdir2diff(node->sd,cent) > maxr2)
306 >                        return;                         /* too far away */
307                  intens = colrval(node->val,GRN) * tri_omegav(node->gdv);
308 <                VSUM(vsum, vsum, node->sd, intens);
308 >                decodedir(sdir, node->sd);
309 >                VSUM(vsum, vsum, sdir, intens);
310                  return;
311          }
312 <        if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
313 <                        DOT(node->gdv[0],cent) > mincos &&
314 <                        DOT(node->gdv[1],cent) > mincos &&
315 <                        DOT(node->gdv[2],cent) > mincos)
316 <                return;
317 <        vector_sum(vsum, &node->kid[0], cent, mincos, ethresh);
318 <        vector_sum(vsum, &node->kid[1], cent, mincos, ethresh);
319 <        vector_sum(vsum, &node->kid[2], cent, mincos, ethresh);
320 <        vector_sum(vsum, &node->kid[3], cent, mincos, ethresh);
312 >        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
313 >                        fdir2diff(node->gdv[0],cent) < maxr2 &&
314 >                        fdir2diff(node->gdv[1],cent) < maxr2 &&
315 >                        fdir2diff(node->gdv[2],cent) < maxr2)
316 >                return;                                 /* containing node */
317 >        vector_sum(vsum, &node->kid[0], cent, maxr2, ethresh);
318 >        vector_sum(vsum, &node->kid[1], cent, maxr2, ethresh);
319 >        vector_sum(vsum, &node->kid[2], cent, maxr2, ethresh);
320 >        vector_sum(vsum, &node->kid[3], cent, maxr2, ethresh);
321   }
322  
323   /* Claim source contributions within the given solid angle */
324   void
325 < claimlight(COLOR intens, TRITREE *node, const FVECT cent, double mincos)
325 > claimlight(COLOR intens, TRITREE *node, FVECT cent, double maxr2)
326   {
327          int     remaining;
328          int     i;
329          if (isleaf(node)) {             /* claim contribution */
330                  COLOR   contrib;
331                  if (node->val[EXP] <= 0)
332 <                        return;
333 <                if (DOT(node->sd,cent) < mincos)
334 <                        return;
332 >                        return;         /* already claimed */
333 >                if (fdir2diff(node->sd,cent) > maxr2)
334 >                        return;         /* too far away */
335                  colr_color(contrib, node->val);
336                  scalecolor(contrib, tri_omegav(node->gdv));
337                  addcolor(intens, contrib);
338                  copycolr(node->val, blkclr);
339                  return;
340          }
341 <        if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
342 <                        DOT(node->gdv[0],cent) > mincos &&
343 <                        DOT(node->gdv[1],cent) > mincos &&
344 <                        DOT(node->gdv[2],cent) > mincos)
345 <                return;
341 >        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
342 >                        fdir2diff(node->gdv[0],cent) < maxr2 &&
343 >                        fdir2diff(node->gdv[1],cent) < maxr2 &&
344 >                        fdir2diff(node->gdv[2],cent) < maxr2)
345 >                return;                 /* previously claimed node */
346          remaining = 0;                  /* recurse on children */
347          for (i = 0; i < 4; i++) {
348 <                claimlight(intens, &node->kid[i], cent, mincos);
348 >                claimlight(intens, &node->kid[i], cent, maxr2);
349                  if (!isleaf(&node->kid[i]) || node->kid[i].val[EXP] != 0)
350                          ++remaining;
351          }
352          if (remaining)
353                  return;
354                                          /* consolidate empties */
355 <        free((void *)node->kid); node->kid = NULL;
355 >        free(node->kid); node->kid = NULL;
356          copycolr(node->val, blkclr);
357 <        VCOPY(node->sd, node->gdv[0]);  /* doesn't really matter */
357 >        node->sd = node->gdv[0];        /* doesn't really matter */
358   }
359  
360   /* Add lost light contribution to the given list */
361   void
362 < add2lost(LOSTLIGHT **llp, COLOR intens, const FVECT cent)
362 > add2lost(LOSTLIGHT **llp, COLOR intens, FVECT cent)
363   {
364          LOSTLIGHT       *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
365  
366          if (newll == NULL)
367                  return;
368          copycolor(newll->intens, intens);
369 <        VCOPY(newll->sd, cent);
369 >        newll->sd = encodedir(cent);
370          newll->next = *llp;
371          *llp = newll;
372   }
373  
374   /* Check lost light list for contributions */
375   void
376 < getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cent, double omega)
376 > getlost(LOSTLIGHT **llp, COLOR intens, FVECT cent, double omega)
377   {
378 <        const double    mincos = 1. - omega/(2.*PI);
378 >        const double    maxr2 = omega/PI;
379          LOSTLIGHT       lhead, *lastp, *thisp;
380  
381          lhead.next = *llp;
382          lastp = &lhead;
383          while ((thisp = lastp->next) != NULL)
384 <                if (DOT(thisp->sd,cent) >= mincos) {
384 >                if (fdir2diff(thisp->sd,cent) <= maxr2) {
385                          LOSTLIGHT       *mynext = thisp->next;
386                          addcolor(intens, thisp->intens);
387 <                        free((void *)thisp);
387 >                        free(thisp);
388                          lastp->next = mynext;
389                  } else
390                          lastp = thisp;
# Line 380 | Line 395 | getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cen
395   void
396   mksources(TRITREE *samptree, double thresh, double maxang)
397   {
398 + #define MAXITER         100
399          const int       ethresh = (int)(log(thresh)/log(2.) + (COLXS+.5));
400          const double    maxomega = 2.*PI*(1. - cos(PI/180./2.*maxang));
401          const double    minintens = .05*thresh*maxomega;
402 +        int             niter = MAXITER;
403          int             nsrcs = 0;
404          LOSTLIGHT       *lostlightlist = NULL;
405          int             emax;
# Line 408 | Line 425 | mksources(TRITREE *samptree, double thresh, double max
425           */
426          if (thresh <= FTINY)
427                  return;
428 <        for ( ; ; ) {
428 >        while (niter--) {
429                  emax = ethresh;         /* find brightest unclaimed */
430                  startleaf = NULL;
431                  for (i = 0; i < NTRUNKBR; i++) {
# Line 419 | Line 436 | mksources(TRITREE *samptree, double thresh, double max
436                  if (startleaf == NULL)
437                          break;
438                                          /* claim it */
439 <                VCOPY(curcent, startleaf->sd);
439 >                decodedir(curcent, startleaf->sd);
440                  curomega = tri_omegav(startleaf->gdv);
441                  currad = sqrt(curomega/PI);
442                  growstep = 3.*currad;
# Line 435 | Line 452 | mksources(TRITREE *samptree, double thresh, double max
452                          vsum[0] = vsum[1] = vsum[2] = .0;
453                          for (i = 0; i < NTRUNKBR; i++)
454                                  vector_sum(vsum, &samptree->kid[i],
455 <                                                curcent, cos(currad+growstep),
455 >                                        curcent, 2.-2.*cos(currad+growstep),
456                                                  thisethresh);
457                          if (normalize(vsum) == .0)
458                                  break;
459 <                        movedist = acos(DOT(vsum,curcent));
459 >                        movedist = Acos(DOT(vsum,curcent));
460                          if (movedist > growstep) {
461                                  VSUB(vsum, vsum, curcent);
462                                  movedist = growstep/VLEN(vsum);
# Line 451 | Line 468 | mksources(TRITREE *samptree, double thresh, double max
468                          curomega = 2.*PI*(1. - cos(currad));
469                          for (i = 0; i < NTRUNKBR; i++)
470                                  claimlight(curintens, &samptree->kid[i],
471 <                                                curcent, cos(currad));
471 >                                                curcent, 2.-2.*cos(currad));
472                  } while (curomega < maxomega &&
473                                  bright(curintens)/curomega > thisthresh);
474                  if (bright(curintens) < minintens) {
# Line 471 | Line 488 | mksources(TRITREE *samptree, double thresh, double max
488                  printf("0\n0\n4 %f %f %f %f\n",
489                                  curcent[0], curcent[1], curcent[2],
490                                  2.*180./PI*currad);
491 +                niter = MAXITER;
492          }
493 + #undef MAXITER
494   }
495  
496   int
# Line 482 | Line 501 | main(int argc, char *argv[])
501          TRITREE *samptree;
502          double  thresh = 0;
503          int     i;
504 <        
505 <        progname = argv[0];
504 >
505 >        fixargv0(argv[0]);              /* sets global progname */
506 >
507          for (i = 1; i < argc && argv[i][0] == '-'; i++)
508                  switch (argv[i][1]) {
509                  case 'd':               /* number of samples */
# Line 536 | Line 556 | userr:
556   }
557  
558   void
559 < eputs(char  *s)
559 > eputs(const char  *s)
560   {
561          static int  midline = 0;
562  
# Line 554 | Line 574 | eputs(char  *s)
574   }
575  
576   void
577 < wputs(char *s)
577 > wputs(const char *s)
578   {
579          /* no warnings */
580   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines