ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mksource.c
Revision: 2.13
Committed: Tue Jun 3 21:31:51 2025 UTC (37 hours, 12 minutes ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.12: +5 -5 lines
Log Message:
refactor: More consistent use of global char * progname and fixargv0()

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mksource.c,v 2.12 2025/01/18 20:23:15 greg Exp $";
3 #endif
4 /*
5 * Generate distant sources corresponding to the given environment map
6 */
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 */
14 #define NTRUNKVERT 4 /* number of vertices at trunk */
15 #define DEF_NSAMPS 262144L /* default # sphere samples */
16 #define DEF_MAXANG 15. /* maximum source angle (deg.) */
17
18 /* Data structure for geodesic samples */
19
20 typedef struct tritree {
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 int32 sd; /* lost source direction */
30 COLOR intens; /* output times solid angle */
31 } LOSTLIGHT;
32
33 FVECT scene_cent; /* center of octree cube */
34 RREAL scene_rad; /* radius to get outside cube from center */
35
36 const COLR blkclr = BLKCOLR;
37
38 #define isleaf(node) ((node)->kid == NULL)
39
40 /* Compute signum of signed volume for three vectors */
41 int
42 vol_sign(const FVECT v1, const FVECT v2, const FVECT v3)
43 {
44 double vol;
45
46 vol = v1[0]*(v2[1]*v3[2] - v2[2]*v3[1]);
47 vol += v1[1]*(v2[2]*v3[0] - v2[0]*v3[2]);
48 vol += v1[2]*(v2[0]*v3[1] - v2[1]*v3[0]);
49 if (vol > .0)
50 return(1);
51 if (vol < .0)
52 return(-1);
53 return(0);
54 }
55
56 /* Is the given direction contained within the specified spherical triangle? */
57 int
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]);
69 if ((sv[0] == sv[1]) & (sv[1] == sv[2]))
70 return(1);
71 return(!sv[0] | !sv[1] | !sv[2]);
72 }
73
74 /* Find leaf containing the given sample direction */
75 TRITREE *
76 findleaf(TRITREE *node, const FVECT sdir)
77 {
78 int i;
79
80 if (isleaf(node))
81 return(intriv(node->gdv,sdir) ? node : (TRITREE *)NULL);
82 for (i = 0; i < 4; i++) {
83 TRITREE *chknode = &node->kid[i];
84 if (intriv(chknode->gdv, sdir))
85 return(isleaf(chknode) ? chknode :
86 findleaf(chknode, sdir));
87 }
88 return(NULL);
89 }
90
91 /* Initialize leaf with random sample inside the given spherical triangle */
92 void
93 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;
101 wt[i] = frandom();
102 j = random() & 1;
103 wt[(i+2-j)%3] = 1. - wt[i] -
104 (wt[(i+1+j)%3] = (1.-wt[i])*frandom());
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, sdir, scene_rad);
115 VCOPY(myray.rdir, sdir);
116 myray.rmax = 0.;
117 ray_trace(&myray);
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, const int32 dv[3])
124 {
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 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], 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 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 */
159 void
160 branchsample(TRITREE *node, int depth)
161 {
162 int i;
163
164 if (depth <= 0)
165 return;
166 if (isleaf(node)) { /* subdivide leaf node */
167 TRITREE branch, *moved_leaf;
168 FVECT sdir;
169 subdivide(&branch, node->gdv);
170 decodedir(sdir, node->sd);
171 moved_leaf = findleaf(&branch, sdir);
172 if (moved_leaf != NULL) { /* bequeath old sample */
173 moved_leaf->sd = node->sd;
174 copycolr(moved_leaf->val, node->val);
175 }
176 for (i = 0; i < 4; i++) /* compute new samples */
177 if (&branch.kid[i] != moved_leaf)
178 leafsample(&branch.kid[i]);
179 *node = branch; /* replace leaf with branch */
180 }
181 for (i = 0; i < 4; i++) /* subdivide children */
182 branchsample(&node->kid[i], depth-1);
183 }
184
185 /* Sample sphere using triangular geodesic mesh */
186 TRITREE *
187 geosample(long nsamps)
188 {
189 int depth;
190 TRITREE *tree;
191 FVECT trunk[NTRUNKVERT];
192 int i, j;
193 /* figure out depth */
194 if ((nsamps -= 4) < 0)
195 error(USER, "minimum number of samples is 4");
196 nsamps = nsamps*(NTRUNKBR-1)/NTRUNKBR; /* round up */
197 for (depth = 0; nsamps > 1; depth++)
198 nsamps >>= 2;
199 /* make base tetrahedron */
200 tree = (TRITREE *)malloc(sizeof(TRITREE));
201 if (tree == NULL) goto memerr;
202 trunk[0][0] = trunk[0][1] = 0; trunk[0][2] = 1;
203 trunk[1][0] = 0;
204 trunk[1][2] = cos(2.*asin(sqrt(2./3.)));
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 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 tree->kid[i].gdv[j] = encodedir(trunk[(i+j)%NTRUNKVERT]);
215 leafsample(&tree->kid[i]);
216 branchsample(&tree->kid[i], depth);
217 }
218 return(tree);
219 memerr:
220 error(SYSTEM, "out of memory in geosample()");
221 return NULL; /* dummy return */
222 }
223
224 /* Compute leaf exponent histogram */
225 void
226 get_ehisto(const TRITREE *node, long exphisto[256])
227 {
228 int i;
229
230 if (isleaf(node)) {
231 ++exphisto[node->val[EXP]];
232 return;
233 }
234 for (i = 0; i < 4; i++)
235 get_ehisto(&node->kid[i], exphisto);
236 }
237
238 /* Get reasonable source threshold */
239 double
240 get_threshold(const TRITREE *tree)
241 {
242 long samptotal = 0;
243 long exphisto[256];
244 int i;
245 /* compute sample histogram */
246 memset(exphisto, 0, sizeof(exphisto));
247 for (i = 0; i < NTRUNKBR; i++)
248 get_ehisto(&tree->kid[i], exphisto);
249 /* use 98th percentile */
250 for (i = 0; i < 256; i++)
251 samptotal += exphisto[i];
252 samptotal /= 50;
253 for (i = 256; (--i > 0) & (samptotal > 0); )
254 samptotal -= exphisto[i];
255 return(ldexp(.75, i-COLXS));
256 }
257
258 /* Find leaf containing the maximum exponent */
259 TRITREE *
260 findemax(TRITREE *node, int *expp)
261 {
262 if (!isleaf(node)) {
263 TRITREE *maxleaf;
264 TRITREE *rleaf;
265 maxleaf = findemax(&node->kid[0], expp);
266 rleaf = findemax(&node->kid[1], expp);
267 if (rleaf != NULL) maxleaf = rleaf;
268 rleaf = findemax(&node->kid[2], expp);
269 if (rleaf != NULL) maxleaf = rleaf;
270 rleaf = findemax(&node->kid[3], expp);
271 if (rleaf != NULL) maxleaf = rleaf;
272 return(maxleaf);
273 }
274 if (node->val[EXP] <= *expp)
275 return(NULL);
276 *expp = node->val[EXP];
277 return(node);
278 }
279
280 /* Compute solid angle of spherical triangle (approx.) */
281 double
282 tri_omegav(const int32 vd[3])
283 {
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 above-threshold perimiter within radius */
296 void
297 vector_sum(FVECT vsum, TRITREE *node,
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; /* 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 decodedir(sdir, node->sd);
309 VSUM(vsum, vsum, sdir, intens);
310 return;
311 }
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, 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; /* 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 (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, 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(node->kid); node->kid = NULL;
356 copycolr(node->val, blkclr);
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, FVECT cent)
363 {
364 LOSTLIGHT *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
365
366 if (newll == NULL)
367 return;
368 copycolor(newll->intens, intens);
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, FVECT cent, double omega)
377 {
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 (fdir2diff(thisp->sd,cent) <= maxr2) {
385 LOSTLIGHT *mynext = thisp->next;
386 addcolor(intens, thisp->intens);
387 free(thisp);
388 lastp->next = mynext;
389 } else
390 lastp = thisp;
391 *llp = lhead.next;
392 }
393
394 /* Create & print distant sources */
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;
406 TRITREE *startleaf;
407 double growstep;
408 FVECT curcent;
409 double currad;
410 double curomega;
411 COLOR curintens;
412 double thisthresh;
413 int thisethresh;
414 int i;
415 /*
416 * General algorithm:
417 * 1) Start with brightest unclaimed pixel
418 * 2) Grow source toward brightest unclaimed perimeter until:
419 * a) Source exceeds maximum size, or
420 * b) Perimeter values all below threshold, or
421 * c) Source average drops below threshold
422 * 3) Loop until nothing over threshold
423 *
424 * Complexity added to absorb insignificant sources in larger ones.
425 */
426 if (thresh <= FTINY)
427 return;
428 while (niter--) {
429 emax = ethresh; /* find brightest unclaimed */
430 startleaf = NULL;
431 for (i = 0; i < NTRUNKBR; i++) {
432 TRITREE *bigger = findemax(&samptree->kid[i], &emax);
433 if (bigger != NULL)
434 startleaf = bigger;
435 }
436 if (startleaf == NULL)
437 break;
438 /* claim it */
439 decodedir(curcent, startleaf->sd);
440 curomega = tri_omegav(startleaf->gdv);
441 currad = sqrt(curomega/PI);
442 growstep = 3.*currad;
443 colr_color(curintens, startleaf->val);
444 thisthresh = .15*bright(curintens);
445 if (thisthresh < thresh) thisthresh = thresh;
446 thisethresh = (int)(log(thisthresh)/log(2.) + (COLXS+.5));
447 scalecolor(curintens, curomega);
448 copycolr(startleaf->val, blkclr);
449 do { /* grow source */
450 FVECT vsum;
451 double movedist;
452 vsum[0] = vsum[1] = vsum[2] = .0;
453 for (i = 0; i < NTRUNKBR; i++)
454 vector_sum(vsum, &samptree->kid[i],
455 curcent, 2.-2.*cos(currad+growstep),
456 thisethresh);
457 if (normalize(vsum) == .0)
458 break;
459 movedist = Acos(DOT(vsum,curcent));
460 if (movedist > growstep) {
461 VSUB(vsum, vsum, curcent);
462 movedist = growstep/VLEN(vsum);
463 VSUM(curcent, curcent, vsum, movedist);
464 normalize(curcent);
465 } else
466 VCOPY(curcent, vsum);
467 currad += growstep;
468 curomega = 2.*PI*(1. - cos(currad));
469 for (i = 0; i < NTRUNKBR; i++)
470 claimlight(curintens, &samptree->kid[i],
471 curcent, 2.-2.*cos(currad));
472 } while (curomega < maxomega &&
473 bright(curintens)/curomega > thisthresh);
474 if (bright(curintens) < minintens) {
475 add2lost(&lostlightlist, curintens, curcent);
476 continue;
477 }
478 /* absorb lost contributions */
479 getlost(&lostlightlist, curintens, curcent, curomega);
480 ++nsrcs; /* output source */
481 scalecolor(curintens, 1./curomega);
482 printf("\nvoid illum IBLout\n");
483 printf("0\n0\n3 %f %f %f\n",
484 colval(curintens,RED),
485 colval(curintens,GRN),
486 colval(curintens,BLU));
487 printf("\nIBLout source IBLsrc%d\n", nsrcs);
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
497 main(int argc, char *argv[])
498 {
499 long nsamps = DEF_NSAMPS;
500 double maxang = DEF_MAXANG;
501 TRITREE *samptree;
502 double thresh = 0;
503 int i;
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 */
510 if (i >= argc-1) goto userr;
511 nsamps = atol(argv[++i]);
512 break;
513 case 't': /* manual threshold */
514 if (i >= argc-1) goto userr;
515 thresh = atof(argv[++i]);
516 break;
517 case 'a': /* maximum source angle */
518 if (i >= argc-1) goto userr;
519 maxang = atof(argv[++i]);
520 if (maxang <= FTINY)
521 goto userr;
522 if (maxang > 180.)
523 maxang = 180.;
524 break;
525 default:
526 goto userr;
527 }
528 if (i < argc-1)
529 goto userr;
530 /* start our ray calculation */
531 directvis = 0;
532 ray_init(i == argc-1 ? argv[i] : (char *)NULL);
533 VCOPY(scene_cent, thescene.cuorg);
534 scene_cent[0] += 0.5*thescene.cusize;
535 scene_cent[1] += 0.5*thescene.cusize;
536 scene_cent[2] += 0.5*thescene.cusize;
537 scene_rad = 0.86603*thescene.cusize;
538 /* sample geodesic mesh */
539 samptree = geosample(nsamps);
540 /* get source threshold */
541 if (thresh <= FTINY)
542 thresh = get_threshold(samptree);
543 /* done with ray samples */
544 ray_done(1);
545 /* print header */
546 printf("# ");
547 printargs(argc, argv, stdout);
548 /* create & print sources */
549 mksources(samptree, thresh, maxang);
550 /* all done, no need to clean up */
551 return(0);
552 userr:
553 fprintf(stderr, "Usage: %s [-d nsamps][-t thresh][-a maxang] [octree]\n",
554 argv[0]);
555 exit(1);
556 }
557
558 void
559 eputs(const char *s)
560 {
561 static int midline = 0;
562
563 if (!*s)
564 return;
565 if (!midline++) {
566 fputs(progname, stderr);
567 fputs(": ", stderr);
568 }
569 fputs(s, stderr);
570 if (s[strlen(s)-1] == '\n') {
571 fflush(stderr);
572 midline = 0;
573 }
574 }
575
576 void
577 wputs(const char *s)
578 {
579 /* no warnings */
580 }