ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mksource.c
Revision: 2.11
Committed: Fri Nov 17 20:02:07 2023 UTC (5 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.10: +2 -4 lines
Log Message:
fix: multiple bug fixes in hyperspectral code, added rvu, mkillum, rsensor, and ranimove to "working except photon map" status

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mksource.c,v 2.10 2023/02/06 22:40:21 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 "resolu.h"
11
12 #define NTRUNKBR 4 /* number of branches at trunk */
13 #define NTRUNKVERT 4 /* number of vertices at trunk */
14 #define DEF_NSAMPS 262144L /* default # sphere samples */
15 #define DEF_MAXANG 15. /* maximum source angle (deg.) */
16
17 /* Data structure for geodesic samples */
18
19 typedef struct tritree {
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 int32 sd; /* lost source direction */
29 COLOR intens; /* output times solid angle */
30 } LOSTLIGHT;
31
32 extern char *progname;
33
34 FVECT scene_cent; /* center of octree cube */
35 RREAL scene_rad; /* radius to get outside cube from center */
36
37 const COLR blkclr = BLKCOLR;
38
39 #define isleaf(node) ((node)->kid == NULL)
40
41 /* Compute signum of signed volume for three vectors */
42 int
43 vol_sign(const FVECT v1, const FVECT v2, const FVECT v3)
44 {
45 double vol;
46
47 vol = v1[0]*(v2[1]*v3[2] - v2[2]*v3[1]);
48 vol += v1[1]*(v2[2]*v3[0] - v2[0]*v3[2]);
49 vol += v1[2]*(v2[0]*v3[1] - v2[1]*v3[0]);
50 if (vol > .0)
51 return(1);
52 if (vol < .0)
53 return(-1);
54 return(0);
55 }
56
57 /* Is the given direction contained within the specified spherical triangle? */
58 int
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]);
70 if ((sv[0] == sv[1]) & (sv[1] == sv[2]))
71 return(1);
72 return(!sv[0] | !sv[1] | !sv[2]);
73 }
74
75 /* Find leaf containing the given sample direction */
76 TRITREE *
77 findleaf(TRITREE *node, const FVECT sdir)
78 {
79 int i;
80
81 if (isleaf(node))
82 return(intriv(node->gdv,sdir) ? node : (TRITREE *)NULL);
83 for (i = 0; i < 4; i++) {
84 TRITREE *chknode = &node->kid[i];
85 if (intriv(chknode->gdv, sdir))
86 return(isleaf(chknode) ? chknode :
87 findleaf(chknode, sdir));
88 }
89 return(NULL);
90 }
91
92 /* Initialize leaf with random sample inside the given spherical triangle */
93 void
94 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;
102 wt[i] = frandom();
103 j = random() & 1;
104 wt[(i+2-j)%3] = 1. - wt[i] -
105 (wt[(i+1+j)%3] = (1.-wt[i])*frandom());
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, sdir, scene_rad);
116 VCOPY(myray.rdir, sdir);
117 myray.rmax = 0.;
118 ray_trace(&myray);
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, const int32 dv[3])
125 {
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 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], 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 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 */
160 void
161 branchsample(TRITREE *node, int depth)
162 {
163 int i;
164
165 if (depth <= 0)
166 return;
167 if (isleaf(node)) { /* subdivide leaf node */
168 TRITREE branch, *moved_leaf;
169 FVECT sdir;
170 subdivide(&branch, node->gdv);
171 decodedir(sdir, node->sd);
172 moved_leaf = findleaf(&branch, sdir);
173 if (moved_leaf != NULL) { /* bequeath old sample */
174 moved_leaf->sd = node->sd;
175 copycolr(moved_leaf->val, node->val);
176 }
177 for (i = 0; i < 4; i++) /* compute new samples */
178 if (&branch.kid[i] != moved_leaf)
179 leafsample(&branch.kid[i]);
180 *node = branch; /* replace leaf with branch */
181 }
182 for (i = 0; i < 4; i++) /* subdivide children */
183 branchsample(&node->kid[i], depth-1);
184 }
185
186 /* Sample sphere using triangular geodesic mesh */
187 TRITREE *
188 geosample(int nsamps)
189 {
190 int depth;
191 TRITREE *tree;
192 FVECT trunk[NTRUNKVERT];
193 int i, j;
194 /* figure out depth */
195 if ((nsamps -= 4) < 0)
196 error(USER, "minimum number of samples is 4");
197 nsamps = nsamps*3/NTRUNKBR; /* round up */
198 for (depth = 0; nsamps > 1; depth++)
199 nsamps >>= 2;
200 /* make base tetrahedron */
201 tree = (TRITREE *)malloc(sizeof(TRITREE));
202 if (tree == NULL) goto memerr;
203 trunk[0][0] = trunk[0][1] = 0; trunk[0][2] = 1;
204 trunk[1][0] = 0;
205 trunk[1][2] = cos(2.*asin(sqrt(2./3.)));
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 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 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 */
226 void
227 get_ehisto(const TRITREE *node, long exphisto[256])
228 {
229 int i;
230
231 if (isleaf(node)) {
232 ++exphisto[node->val[EXP]];
233 return;
234 }
235 for (i = 0; i < 4; i++)
236 get_ehisto(&node->kid[i], exphisto);
237 }
238
239 /* Get reasonable source threshold */
240 double
241 get_threshold(const TRITREE *tree)
242 {
243 long samptotal = 0;
244 long exphisto[256];
245 int i;
246 /* compute sample histogram */
247 memset(exphisto, 0, sizeof(exphisto));
248 for (i = 0; i < NTRUNKBR; i++)
249 get_ehisto(&tree->kid[i], exphisto);
250 /* use 98th percentile */
251 for (i = 0; i < 256; i++)
252 samptotal += exphisto[i];
253 samptotal /= 50;
254 for (i = 256; (--i > 0) & (samptotal > 0); )
255 samptotal -= exphisto[i];
256 return(ldexp(.75, i-COLXS));
257 }
258
259 /* Find leaf containing the maximum exponent */
260 TRITREE *
261 findemax(TRITREE *node, int *expp)
262 {
263 if (!isleaf(node)) {
264 TRITREE *maxleaf;
265 TRITREE *rleaf;
266 maxleaf = findemax(&node->kid[0], expp);
267 rleaf = findemax(&node->kid[1], expp);
268 if (rleaf != NULL) maxleaf = rleaf;
269 rleaf = findemax(&node->kid[2], expp);
270 if (rleaf != NULL) maxleaf = rleaf;
271 rleaf = findemax(&node->kid[3], expp);
272 if (rleaf != NULL) maxleaf = rleaf;
273 return(maxleaf);
274 }
275 if (node->val[EXP] <= *expp)
276 return(NULL);
277 *expp = node->val[EXP];
278 return(node);
279 }
280
281 /* Compute solid angle of spherical triangle (approx.) */
282 double
283 tri_omegav(const int32 vd[3])
284 {
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 above-threshold perimiter within radius */
297 void
298 vector_sum(FVECT vsum, TRITREE *node,
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; /* 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 decodedir(sdir, node->sd);
310 VSUM(vsum, vsum, sdir, intens);
311 return;
312 }
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, 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; /* 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 (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, 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(node->kid); node->kid = NULL;
357 copycolr(node->val, blkclr);
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, FVECT cent)
364 {
365 LOSTLIGHT *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
366
367 if (newll == NULL)
368 return;
369 copycolor(newll->intens, intens);
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, FVECT cent, double omega)
378 {
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 (fdir2diff(thisp->sd,cent) <= maxr2) {
386 LOSTLIGHT *mynext = thisp->next;
387 addcolor(intens, thisp->intens);
388 free(thisp);
389 lastp->next = mynext;
390 } else
391 lastp = thisp;
392 *llp = lhead.next;
393 }
394
395 /* Create & print distant sources */
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;
408 double growstep;
409 FVECT curcent;
410 double currad;
411 double curomega;
412 COLOR curintens;
413 double thisthresh;
414 int thisethresh;
415 int i;
416 /*
417 * General algorithm:
418 * 1) Start with brightest unclaimed pixel
419 * 2) Grow source toward brightest unclaimed perimeter until:
420 * a) Source exceeds maximum size, or
421 * b) Perimeter values all below threshold, or
422 * c) Source average drops below threshold
423 * 3) Loop until nothing over threshold
424 *
425 * Complexity added to absorb insignificant sources in larger ones.
426 */
427 if (thresh <= FTINY)
428 return;
429 while (niter--) {
430 emax = ethresh; /* find brightest unclaimed */
431 startleaf = NULL;
432 for (i = 0; i < NTRUNKBR; i++) {
433 TRITREE *bigger = findemax(&samptree->kid[i], &emax);
434 if (bigger != NULL)
435 startleaf = bigger;
436 }
437 if (startleaf == NULL)
438 break;
439 /* claim it */
440 decodedir(curcent, startleaf->sd);
441 curomega = tri_omegav(startleaf->gdv);
442 currad = sqrt(curomega/PI);
443 growstep = 3.*currad;
444 colr_color(curintens, startleaf->val);
445 thisthresh = .15*bright(curintens);
446 if (thisthresh < thresh) thisthresh = thresh;
447 thisethresh = (int)(log(thisthresh)/log(2.) + (COLXS+.5));
448 scalecolor(curintens, curomega);
449 copycolr(startleaf->val, blkclr);
450 do { /* grow source */
451 FVECT vsum;
452 double movedist;
453 vsum[0] = vsum[1] = vsum[2] = .0;
454 for (i = 0; i < NTRUNKBR; i++)
455 vector_sum(vsum, &samptree->kid[i],
456 curcent, 2.-2.*cos(currad+growstep),
457 thisethresh);
458 if (normalize(vsum) == .0)
459 break;
460 movedist = Acos(DOT(vsum,curcent));
461 if (movedist > growstep) {
462 VSUB(vsum, vsum, curcent);
463 movedist = growstep/VLEN(vsum);
464 VSUM(curcent, curcent, vsum, movedist);
465 normalize(curcent);
466 } else
467 VCOPY(curcent, vsum);
468 currad += growstep;
469 curomega = 2.*PI*(1. - cos(currad));
470 for (i = 0; i < NTRUNKBR; i++)
471 claimlight(curintens, &samptree->kid[i],
472 curcent, 2.-2.*cos(currad));
473 } while (curomega < maxomega &&
474 bright(curintens)/curomega > thisthresh);
475 if (bright(curintens) < minintens) {
476 add2lost(&lostlightlist, curintens, curcent);
477 continue;
478 }
479 /* absorb lost contributions */
480 getlost(&lostlightlist, curintens, curcent, curomega);
481 ++nsrcs; /* output source */
482 scalecolor(curintens, 1./curomega);
483 printf("\nvoid illum IBLout\n");
484 printf("0\n0\n3 %f %f %f\n",
485 colval(curintens,RED),
486 colval(curintens,GRN),
487 colval(curintens,BLU));
488 printf("\nIBLout source IBLsrc%d\n", nsrcs);
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
498 main(int argc, char *argv[])
499 {
500 long nsamps = DEF_NSAMPS;
501 double maxang = DEF_MAXANG;
502 TRITREE *samptree;
503 double thresh = 0;
504 int i;
505
506 progname = argv[0];
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 }