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.2 by greg, Sat Jul 30 17:01:30 2005 UTC vs.
Revision 2.11 by greg, Fri Nov 17 20:02:07 2023 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7  
8   #include "ray.h"
9   #include "random.h"
10 + #include "resolu.h"
11  
12   #define NTRUNKBR        4               /* number of branches at trunk */
13   #define NTRUNKVERT      4               /* number of vertices at trunk */
# Line 16 | 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  
32 < char    *progname;
32 > extern char     *progname;
33  
34   FVECT   scene_cent;             /* center of octree cube */
35   RREAL   scene_rad;              /* radius to get outside cube from center */
# Line 55 | 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 90 | 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 97 | 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),
110 <                        colval(myray.rcol,GRN),
111 <                        colval(myray.rcol,BLU));
119 >        scolor_colr(leaf->val, myray.rcol);
120   }
121  
122   /* Initialize a branch node contained in the given spherical triangle */
123   void
124 < subdivide(TRITREE *branch, FVECT dv[3])
124 > subdivide(TRITREE *branch, const int32 dv[3])
125   {
126 <        FVECT   sdv[3];
126 >        FVECT   dvv[3], sdv[3];
127 >        int32   sd[3];
128          int     i;
129          
130 <        for (i = 0; i < 3; i++)         /* copy spherical triangle */
131 <                VCOPY(branch->gdv[i], dv[i]);
130 >        for (i = 0; i < 3; i++) {       /* copy spherical triangle */
131 >                branch->gdv[i] = dv[i];
132 >                decodedir(dvv[i], dv[i]);
133 >        }
134          for (i = 0; i < 3; i++) {       /* create new vertices */
135                  int     j = (i+1)%3;
136 <                VADD(sdv[i], dv[i], dv[j]);
136 >                VADD(sdv[i], dvv[i], dvv[j]);
137                  normalize(sdv[i]);
138 +                sd[i] = encodedir(sdv[i]);
139          }
140                                          /* allocate leaves */
141          branch->kid = (TRITREE *)calloc(4, sizeof(TRITREE));
142          if (branch->kid == NULL)
143                  error(SYSTEM, "out of memory in subdivide()");
144                                          /* assign subtriangle directions */
145 <        VCOPY(branch->kid[0].gdv[0], dv[0]);
146 <        VCOPY(branch->kid[0].gdv[1], sdv[0]);
147 <        VCOPY(branch->kid[0].gdv[2], sdv[2]);
148 <        VCOPY(branch->kid[1].gdv[0], sdv[0]);
149 <        VCOPY(branch->kid[1].gdv[1], dv[1]);
150 <        VCOPY(branch->kid[1].gdv[2], sdv[1]);
151 <        VCOPY(branch->kid[2].gdv[0], sdv[1]);
152 <        VCOPY(branch->kid[2].gdv[1], dv[2]);
153 <        VCOPY(branch->kid[2].gdv[2], sdv[2]);
154 <        VCOPY(branch->kid[3].gdv[0], sdv[0]);
155 <        VCOPY(branch->kid[3].gdv[1], sdv[1]);
156 <        VCOPY(branch->kid[3].gdv[2], sdv[2]);
145 >        branch->kid[0].gdv[0] = dv[0];
146 >        branch->kid[0].gdv[1] = sd[0];
147 >        branch->kid[0].gdv[2] = sd[2];
148 >        branch->kid[1].gdv[0] = sd[0];
149 >        branch->kid[1].gdv[1] = dv[1];
150 >        branch->kid[1].gdv[2] = sd[1];
151 >        branch->kid[2].gdv[0] = sd[1];
152 >        branch->kid[2].gdv[1] = dv[2];
153 >        branch->kid[2].gdv[2] = sd[2];
154 >        branch->kid[3].gdv[0] = sd[0];
155 >        branch->kid[3].gdv[1] = sd[1];
156 >        branch->kid[3].gdv[2] = sd[2];
157   }
158  
159   /* Recursively subdivide the given node to the specified quadtree depth */
# Line 154 | Line 166 | branchsample(TRITREE *node, int depth)
166                  return;
167          if (isleaf(node)) {                     /* subdivide leaf node */
168                  TRITREE branch, *moved_leaf;
169 +                FVECT   sdir;
170                  subdivide(&branch, node->gdv);
171 <                moved_leaf = findleaf(&branch, node->sd);
171 >                decodedir(sdir, node->sd);
172 >                moved_leaf = findleaf(&branch, sdir);
173                  if (moved_leaf != NULL) {       /* bequeath old sample */
174 <                        VCOPY(moved_leaf->sd, node->sd);
174 >                        moved_leaf->sd = node->sd;
175                          copycolr(moved_leaf->val, node->val);
176                  }
177                  for (i = 0; i < 4; i++)         /* compute new samples */
# Line 192 | Line 206 | geosample(int nsamps)
206          trunk[1][1] = sqrt(1. - trunk[1][2]*trunk[1][2]);
207          spinvector(trunk[2], trunk[1], trunk[0], 2.*PI/3.);
208          spinvector(trunk[3], trunk[1], trunk[0], 4.*PI/3.);
209 <        VCOPY(tree->gdv[0], trunk[0]);
196 <        VCOPY(tree->gdv[1], trunk[0]);
197 <        VCOPY(tree->gdv[2], trunk[0]);
209 >        tree->gdv[0] = tree->gdv[1] = tree->gdv[2] = encodedir(trunk[0]);
210          tree->kid = (TRITREE *)calloc(NTRUNKBR, sizeof(TRITREE));
211          if (tree->kid == NULL) goto memerr;
212                                          /* grow our tree from trunk */
213          for (i = 0; i < NTRUNKBR; i++) {
214                  for (j = 0; j < 3; j++)         /* XXX works for tetra only */
215 <                        VCOPY(tree->kid[i].gdv[j], trunk[(i+j)%NTRUNKVERT]);
215 >                        tree->kid[i].gdv[j] = encodedir(trunk[(i+j)%NTRUNKVERT]);
216                  leafsample(&tree->kid[i]);
217                  branchsample(&tree->kid[i], depth);
218          }
219          return(tree);
220   memerr:
221          error(SYSTEM, "out of memory in geosample()");
222 +        return NULL; /* dummy return */
223   }
224  
225   /* Compute leaf exponent histogram */
# Line 227 | Line 240 | get_ehisto(const TRITREE *node, long exphisto[256])
240   double
241   get_threshold(const TRITREE *tree)
242   {
243 +        long    samptotal = 0;
244          long    exphisto[256];
231        long    samptotal;
245          int     i;
246                                                  /* compute sample histogram */
247 <        memset((void *)exphisto, 0, sizeof(exphisto));
247 >        memset(exphisto, 0, sizeof(exphisto));
248          for (i = 0; i < NTRUNKBR; i++)
249                  get_ehisto(&tree->kid[i], exphisto);
250                                                  /* use 98th percentile */
# Line 267 | Line 280 | findemax(TRITREE *node, int *expp)
280  
281   /* Compute solid angle of spherical triangle (approx.) */
282   double
283 < tri_omegav(FVECT v[3])
283 > tri_omegav(const int32 vd[3])
284   {
285 <        FVECT   e1, e2, vcross;
286 <        
285 >        FVECT   v[3], e1, e2, vcross;
286 >
287 >        decodedir(v[0], vd[0]);
288 >        decodedir(v[1], vd[1]);
289 >        decodedir(v[2], vd[2]);
290          VSUB(e1, v[1], v[0]);
291          VSUB(e2, v[2], v[1]);
292          fcross(vcross, e1, e2);
293          return(.5*VLEN(vcross));
294   }
295  
296 < /* Sum intensity times direction for non-zero leaves */
296 > /* Sum intensity times direction for above-threshold perimiter within radius */
297   void
298   vector_sum(FVECT vsum, TRITREE *node,
299 <                const FVECT cent, double mincos, int ethresh)
299 >                FVECT cent, double maxr2, int ethresh)
300   {
301          if (isleaf(node)) {
302                  double  intens;
303 +                FVECT   sdir;
304                  if (node->val[EXP] < ethresh)
305 <                        return;
306 <                if (DOT(node->sd,cent) < mincos)
307 <                        return;
305 >                        return;                         /* below threshold */
306 >                if (fdir2diff(node->sd,cent) > maxr2)
307 >                        return;                         /* too far away */
308                  intens = colrval(node->val,GRN) * tri_omegav(node->gdv);
309 <                VSUM(vsum, vsum, node->sd, intens);
309 >                decodedir(sdir, node->sd);
310 >                VSUM(vsum, vsum, sdir, intens);
311                  return;
312          }
313 <        if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
314 <                        DOT(node->gdv[0],cent) > mincos &&
315 <                        DOT(node->gdv[1],cent) > mincos &&
316 <                        DOT(node->gdv[2],cent) > mincos)
317 <                return;
318 <        vector_sum(vsum, &node->kid[0], cent, mincos, ethresh);
319 <        vector_sum(vsum, &node->kid[1], cent, mincos, ethresh);
320 <        vector_sum(vsum, &node->kid[2], cent, mincos, ethresh);
321 <        vector_sum(vsum, &node->kid[3], cent, mincos, ethresh);
313 >        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
314 >                        fdir2diff(node->gdv[0],cent) < maxr2 &&
315 >                        fdir2diff(node->gdv[1],cent) < maxr2 &&
316 >                        fdir2diff(node->gdv[2],cent) < maxr2)
317 >                return;                                 /* containing node */
318 >        vector_sum(vsum, &node->kid[0], cent, maxr2, ethresh);
319 >        vector_sum(vsum, &node->kid[1], cent, maxr2, ethresh);
320 >        vector_sum(vsum, &node->kid[2], cent, maxr2, ethresh);
321 >        vector_sum(vsum, &node->kid[3], cent, maxr2, ethresh);
322   }
323  
324   /* Claim source contributions within the given solid angle */
325   void
326 < claimlight(COLOR intens, TRITREE *node, const FVECT cent, double mincos)
326 > claimlight(COLOR intens, TRITREE *node, FVECT cent, double maxr2)
327   {
328          int     remaining;
329          int     i;
330          if (isleaf(node)) {             /* claim contribution */
331                  COLOR   contrib;
332                  if (node->val[EXP] <= 0)
333 <                        return;
334 <                if (DOT(node->sd,cent) < mincos)
335 <                        return;
333 >                        return;         /* already claimed */
334 >                if (fdir2diff(node->sd,cent) > maxr2)
335 >                        return;         /* too far away */
336                  colr_color(contrib, node->val);
337                  scalecolor(contrib, tri_omegav(node->gdv));
338                  addcolor(intens, contrib);
339                  copycolr(node->val, blkclr);
340                  return;
341          }
342 <        if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
343 <                        DOT(node->gdv[0],cent) > mincos &&
344 <                        DOT(node->gdv[1],cent) > mincos &&
345 <                        DOT(node->gdv[2],cent) > mincos)
346 <                return;
342 >        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
343 >                        fdir2diff(node->gdv[0],cent) < maxr2 &&
344 >                        fdir2diff(node->gdv[1],cent) < maxr2 &&
345 >                        fdir2diff(node->gdv[2],cent) < maxr2)
346 >                return;                 /* previously claimed node */
347          remaining = 0;                  /* recurse on children */
348          for (i = 0; i < 4; i++) {
349 <                claimlight(intens, &node->kid[i], cent, mincos);
349 >                claimlight(intens, &node->kid[i], cent, maxr2);
350                  if (!isleaf(&node->kid[i]) || node->kid[i].val[EXP] != 0)
351                          ++remaining;
352          }
353          if (remaining)
354                  return;
355                                          /* consolidate empties */
356 <        free((void *)node->kid); node->kid = NULL;
356 >        free(node->kid); node->kid = NULL;
357          copycolr(node->val, blkclr);
358 <        VCOPY(node->sd, node->gdv[0]);  /* doesn't really matter */
358 >        node->sd = node->gdv[0];        /* doesn't really matter */
359   }
360  
361   /* Add lost light contribution to the given list */
362   void
363 < add2lost(LOSTLIGHT **llp, COLOR intens, const FVECT cent)
363 > add2lost(LOSTLIGHT **llp, COLOR intens, FVECT cent)
364   {
365          LOSTLIGHT       *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
366  
367          if (newll == NULL)
368                  return;
369          copycolor(newll->intens, intens);
370 <        VCOPY(newll->sd, cent);
370 >        newll->sd = encodedir(cent);
371          newll->next = *llp;
372          *llp = newll;
373   }
374  
375   /* Check lost light list for contributions */
376   void
377 < getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cent, double omega)
377 > getlost(LOSTLIGHT **llp, COLOR intens, FVECT cent, double omega)
378   {
379 <        const double    mincos = 1. - omega/(2.*PI);
379 >        const double    maxr2 = omega/PI;
380          LOSTLIGHT       lhead, *lastp, *thisp;
381  
382          lhead.next = *llp;
383          lastp = &lhead;
384          while ((thisp = lastp->next) != NULL)
385 <                if (DOT(thisp->sd,cent) >= mincos) {
385 >                if (fdir2diff(thisp->sd,cent) <= maxr2) {
386                          LOSTLIGHT       *mynext = thisp->next;
387                          addcolor(intens, thisp->intens);
388 <                        free((void *)thisp);
388 >                        free(thisp);
389                          lastp->next = mynext;
390                  } else
391                          lastp = thisp;
# Line 378 | Line 396 | getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cen
396   void
397   mksources(TRITREE *samptree, double thresh, double maxang)
398   {
399 + #define MAXITER         100
400          const int       ethresh = (int)(log(thresh)/log(2.) + (COLXS+.5));
401          const double    maxomega = 2.*PI*(1. - cos(PI/180./2.*maxang));
402          const double    minintens = .05*thresh*maxomega;
403 +        int             niter = MAXITER;
404          int             nsrcs = 0;
405          LOSTLIGHT       *lostlightlist = NULL;
406          int             emax;
407          TRITREE         *startleaf;
388        COLOR           cval;
408          double          growstep;
409          FVECT           curcent;
410          double          currad;
# Line 407 | Line 426 | mksources(TRITREE *samptree, double thresh, double max
426           */
427          if (thresh <= FTINY)
428                  return;
429 <        for ( ; ; ) {
429 >        while (niter--) {
430                  emax = ethresh;         /* find brightest unclaimed */
431                  startleaf = NULL;
432                  for (i = 0; i < NTRUNKBR; i++) {
# Line 418 | 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 434 | 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;
460 <                        movedist = acos(DOT(vsum,curcent));
460 >                        movedist = Acos(DOT(vsum,curcent));
461                          if (movedist > growstep) {
462                                  VSUB(vsum, vsum, curcent);
463                                  movedist = growstep/VLEN(vsum);
# Line 450 | 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) {
# Line 470 | Line 489 | mksources(TRITREE *samptree, double thresh, double max
489                  printf("0\n0\n4 %f %f %f %f\n",
490                                  curcent[0], curcent[1], curcent[2],
491                                  2.*180./PI*currad);
492 +                niter = MAXITER;
493          }
494 + #undef MAXITER
495   }
496  
497   int
# Line 535 | Line 556 | userr:
556   }
557  
558   void
559 < eputs(char  *s)
559 > eputs(const char  *s)
560   {
561          static int  midline = 0;
562  
# Line 553 | 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