ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mksource.c
Revision: 2.2
Committed: Sat Jul 30 17:01:30 2005 UTC (18 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1
Changes since 2.1: +3 -3 lines
Log Message:
Changed to 98th percentile for default threshold

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mksource.c,v 2.1 2005/04/12 03:30:43 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
11 #define NTRUNKBR 4 /* number of branches at trunk */
12 #define NTRUNKVERT 4 /* number of vertices at trunk */
13 #define DEF_NSAMPS 262144L /* default # sphere samples */
14 #define DEF_MAXANG 15. /* maximum source angle (deg.) */
15
16 /* Data structure for geodesic samples */
17
18 typedef struct tritree {
19 FVECT gdv[3]; /* spherical triangle vertex direc. */
20 FVECT sd; /* sample direction if leaf */
21 struct tritree *kid; /* 4 children if branch node */
22 COLR val; /* sampled color value */
23 } TRITREE;
24
25 typedef struct lostlight {
26 struct lostlight *next; /* next in list */
27 FVECT sd; /* lost source direction */
28 COLOR intens; /* output times solid angle */
29 } LOSTLIGHT;
30
31 char *progname;
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(FVECT tri[3], const FVECT sdir)
59 {
60 int sv[3];
61
62 sv[0] = vol_sign(sdir, tri[0], tri[1]);
63 sv[1] = vol_sign(sdir, tri[1], tri[2]);
64 sv[2] = vol_sign(sdir, tri[2], tri[0]);
65 if ((sv[0] == sv[1]) & (sv[1] == sv[2]))
66 return(1);
67 return(!sv[0] | !sv[1] | !sv[2]);
68 }
69
70 /* Find leaf containing the given sample direction */
71 TRITREE *
72 findleaf(TRITREE *node, const FVECT sdir)
73 {
74 int i;
75
76 if (isleaf(node))
77 return(intriv(node->gdv,sdir) ? node : (TRITREE *)NULL);
78 for (i = 0; i < 4; i++) {
79 TRITREE *chknode = &node->kid[i];
80 if (intriv(chknode->gdv, sdir))
81 return(isleaf(chknode) ? chknode :
82 findleaf(chknode, sdir));
83 }
84 return(NULL);
85 }
86
87 /* Initialize leaf with random sample inside the given spherical triangle */
88 void
89 leafsample(TRITREE *leaf)
90 {
91 RAY myray;
92 RREAL wt[3];
93 int i, j;
94 /* random point on triangle */
95 i = random() % 3;
96 wt[i] = frandom();
97 j = random() & 1;
98 wt[(i+2-j)%3] = 1. - wt[i] -
99 (wt[(i+1+j)%3] = (1.-wt[i])*frandom());
100 leaf->sd[0] = leaf->sd[1] = leaf->sd[2] = .0;
101 for (i = 0; i < 3; i++)
102 VSUM(leaf->sd, leaf->sd, leaf->gdv[i], wt[i]);
103 normalize(leaf->sd); /* record sample direction */
104 /* evaluate at inf. */
105 VSUM(myray.rorg, scene_cent, leaf->sd, scene_rad);
106 VCOPY(myray.rdir, leaf->sd);
107 myray.rmax = 0.;
108 ray_trace(&myray);
109 setcolr(leaf->val, colval(myray.rcol,RED),
110 colval(myray.rcol,GRN),
111 colval(myray.rcol,BLU));
112 }
113
114 /* Initialize a branch node contained in the given spherical triangle */
115 void
116 subdivide(TRITREE *branch, FVECT dv[3])
117 {
118 FVECT sdv[3];
119 int i;
120
121 for (i = 0; i < 3; i++) /* copy spherical triangle */
122 VCOPY(branch->gdv[i], dv[i]);
123 for (i = 0; i < 3; i++) { /* create new vertices */
124 int j = (i+1)%3;
125 VADD(sdv[i], dv[i], dv[j]);
126 normalize(sdv[i]);
127 }
128 /* allocate leaves */
129 branch->kid = (TRITREE *)calloc(4, sizeof(TRITREE));
130 if (branch->kid == NULL)
131 error(SYSTEM, "out of memory in subdivide()");
132 /* assign subtriangle directions */
133 VCOPY(branch->kid[0].gdv[0], dv[0]);
134 VCOPY(branch->kid[0].gdv[1], sdv[0]);
135 VCOPY(branch->kid[0].gdv[2], sdv[2]);
136 VCOPY(branch->kid[1].gdv[0], sdv[0]);
137 VCOPY(branch->kid[1].gdv[1], dv[1]);
138 VCOPY(branch->kid[1].gdv[2], sdv[1]);
139 VCOPY(branch->kid[2].gdv[0], sdv[1]);
140 VCOPY(branch->kid[2].gdv[1], dv[2]);
141 VCOPY(branch->kid[2].gdv[2], sdv[2]);
142 VCOPY(branch->kid[3].gdv[0], sdv[0]);
143 VCOPY(branch->kid[3].gdv[1], sdv[1]);
144 VCOPY(branch->kid[3].gdv[2], sdv[2]);
145 }
146
147 /* Recursively subdivide the given node to the specified quadtree depth */
148 void
149 branchsample(TRITREE *node, int depth)
150 {
151 int i;
152
153 if (depth <= 0)
154 return;
155 if (isleaf(node)) { /* subdivide leaf node */
156 TRITREE branch, *moved_leaf;
157 subdivide(&branch, node->gdv);
158 moved_leaf = findleaf(&branch, node->sd);
159 if (moved_leaf != NULL) { /* bequeath old sample */
160 VCOPY(moved_leaf->sd, node->sd);
161 copycolr(moved_leaf->val, node->val);
162 }
163 for (i = 0; i < 4; i++) /* compute new samples */
164 if (&branch.kid[i] != moved_leaf)
165 leafsample(&branch.kid[i]);
166 *node = branch; /* replace leaf with branch */
167 }
168 for (i = 0; i < 4; i++) /* subdivide children */
169 branchsample(&node->kid[i], depth-1);
170 }
171
172 /* Sample sphere using triangular geodesic mesh */
173 TRITREE *
174 geosample(int nsamps)
175 {
176 int depth;
177 TRITREE *tree;
178 FVECT trunk[NTRUNKVERT];
179 int i, j;
180 /* figure out depth */
181 if ((nsamps -= 4) < 0)
182 error(USER, "minimum number of samples is 4");
183 nsamps = nsamps*3/NTRUNKBR; /* round up */
184 for (depth = 0; nsamps > 1; depth++)
185 nsamps >>= 2;
186 /* make base tetrahedron */
187 tree = (TRITREE *)malloc(sizeof(TRITREE));
188 if (tree == NULL) goto memerr;
189 trunk[0][0] = trunk[0][1] = 0; trunk[0][2] = 1;
190 trunk[1][0] = 0;
191 trunk[1][2] = cos(2.*asin(sqrt(2./3.)));
192 trunk[1][1] = sqrt(1. - trunk[1][2]*trunk[1][2]);
193 spinvector(trunk[2], trunk[1], trunk[0], 2.*PI/3.);
194 spinvector(trunk[3], trunk[1], trunk[0], 4.*PI/3.);
195 VCOPY(tree->gdv[0], trunk[0]);
196 VCOPY(tree->gdv[1], trunk[0]);
197 VCOPY(tree->gdv[2], trunk[0]);
198 tree->kid = (TRITREE *)calloc(NTRUNKBR, sizeof(TRITREE));
199 if (tree->kid == NULL) goto memerr;
200 /* grow our tree from trunk */
201 for (i = 0; i < NTRUNKBR; i++) {
202 for (j = 0; j < 3; j++) /* XXX works for tetra only */
203 VCOPY(tree->kid[i].gdv[j], trunk[(i+j)%NTRUNKVERT]);
204 leafsample(&tree->kid[i]);
205 branchsample(&tree->kid[i], depth);
206 }
207 return(tree);
208 memerr:
209 error(SYSTEM, "out of memory in geosample()");
210 }
211
212 /* Compute leaf exponent histogram */
213 void
214 get_ehisto(const TRITREE *node, long exphisto[256])
215 {
216 int i;
217
218 if (isleaf(node)) {
219 ++exphisto[node->val[EXP]];
220 return;
221 }
222 for (i = 0; i < 4; i++)
223 get_ehisto(&node->kid[i], exphisto);
224 }
225
226 /* Get reasonable source threshold */
227 double
228 get_threshold(const TRITREE *tree)
229 {
230 long exphisto[256];
231 long samptotal;
232 int i;
233 /* compute sample histogram */
234 memset((void *)exphisto, 0, sizeof(exphisto));
235 for (i = 0; i < NTRUNKBR; i++)
236 get_ehisto(&tree->kid[i], exphisto);
237 /* use 98th percentile */
238 for (i = 0; i < 256; i++)
239 samptotal += exphisto[i];
240 samptotal /= 50;
241 for (i = 256; (--i > 0) & (samptotal > 0); )
242 samptotal -= exphisto[i];
243 return(ldexp(.75, i-COLXS));
244 }
245
246 /* Find leaf containing the maximum exponent */
247 TRITREE *
248 findemax(TRITREE *node, int *expp)
249 {
250 if (!isleaf(node)) {
251 TRITREE *maxleaf;
252 TRITREE *rleaf;
253 maxleaf = findemax(&node->kid[0], expp);
254 rleaf = findemax(&node->kid[1], expp);
255 if (rleaf != NULL) maxleaf = rleaf;
256 rleaf = findemax(&node->kid[2], expp);
257 if (rleaf != NULL) maxleaf = rleaf;
258 rleaf = findemax(&node->kid[3], expp);
259 if (rleaf != NULL) maxleaf = rleaf;
260 return(maxleaf);
261 }
262 if (node->val[EXP] <= *expp)
263 return(NULL);
264 *expp = node->val[EXP];
265 return(node);
266 }
267
268 /* Compute solid angle of spherical triangle (approx.) */
269 double
270 tri_omegav(FVECT v[3])
271 {
272 FVECT e1, e2, vcross;
273
274 VSUB(e1, v[1], v[0]);
275 VSUB(e2, v[2], v[1]);
276 fcross(vcross, e1, e2);
277 return(.5*VLEN(vcross));
278 }
279
280 /* Sum intensity times direction for non-zero leaves */
281 void
282 vector_sum(FVECT vsum, TRITREE *node,
283 const FVECT cent, double mincos, int ethresh)
284 {
285 if (isleaf(node)) {
286 double intens;
287 if (node->val[EXP] < ethresh)
288 return;
289 if (DOT(node->sd,cent) < mincos)
290 return;
291 intens = colrval(node->val,GRN) * tri_omegav(node->gdv);
292 VSUM(vsum, vsum, node->sd, intens);
293 return;
294 }
295 if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
296 DOT(node->gdv[0],cent) > mincos &&
297 DOT(node->gdv[1],cent) > mincos &&
298 DOT(node->gdv[2],cent) > mincos)
299 return;
300 vector_sum(vsum, &node->kid[0], cent, mincos, ethresh);
301 vector_sum(vsum, &node->kid[1], cent, mincos, ethresh);
302 vector_sum(vsum, &node->kid[2], cent, mincos, ethresh);
303 vector_sum(vsum, &node->kid[3], cent, mincos, ethresh);
304 }
305
306 /* Claim source contributions within the given solid angle */
307 void
308 claimlight(COLOR intens, TRITREE *node, const FVECT cent, double mincos)
309 {
310 int remaining;
311 int i;
312 if (isleaf(node)) { /* claim contribution */
313 COLOR contrib;
314 if (node->val[EXP] <= 0)
315 return;
316 if (DOT(node->sd,cent) < mincos)
317 return;
318 colr_color(contrib, node->val);
319 scalecolor(contrib, tri_omegav(node->gdv));
320 addcolor(intens, contrib);
321 copycolr(node->val, blkclr);
322 return;
323 }
324 if (DOT(node->gdv[0],node->gdv[1]) < mincos &&
325 DOT(node->gdv[0],cent) > mincos &&
326 DOT(node->gdv[1],cent) > mincos &&
327 DOT(node->gdv[2],cent) > mincos)
328 return;
329 remaining = 0; /* recurse on children */
330 for (i = 0; i < 4; i++) {
331 claimlight(intens, &node->kid[i], cent, mincos);
332 if (!isleaf(&node->kid[i]) || node->kid[i].val[EXP] != 0)
333 ++remaining;
334 }
335 if (remaining)
336 return;
337 /* consolidate empties */
338 free((void *)node->kid); node->kid = NULL;
339 copycolr(node->val, blkclr);
340 VCOPY(node->sd, node->gdv[0]); /* doesn't really matter */
341 }
342
343 /* Add lost light contribution to the given list */
344 void
345 add2lost(LOSTLIGHT **llp, COLOR intens, const FVECT cent)
346 {
347 LOSTLIGHT *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
348
349 if (newll == NULL)
350 return;
351 copycolor(newll->intens, intens);
352 VCOPY(newll->sd, cent);
353 newll->next = *llp;
354 *llp = newll;
355 }
356
357 /* Check lost light list for contributions */
358 void
359 getlost(LOSTLIGHT **llp, COLOR intens, const FVECT cent, double omega)
360 {
361 const double mincos = 1. - omega/(2.*PI);
362 LOSTLIGHT lhead, *lastp, *thisp;
363
364 lhead.next = *llp;
365 lastp = &lhead;
366 while ((thisp = lastp->next) != NULL)
367 if (DOT(thisp->sd,cent) >= mincos) {
368 LOSTLIGHT *mynext = thisp->next;
369 addcolor(intens, thisp->intens);
370 free((void *)thisp);
371 lastp->next = mynext;
372 } else
373 lastp = thisp;
374 *llp = lhead.next;
375 }
376
377 /* Create & print distant sources */
378 void
379 mksources(TRITREE *samptree, double thresh, double maxang)
380 {
381 const int ethresh = (int)(log(thresh)/log(2.) + (COLXS+.5));
382 const double maxomega = 2.*PI*(1. - cos(PI/180./2.*maxang));
383 const double minintens = .05*thresh*maxomega;
384 int nsrcs = 0;
385 LOSTLIGHT *lostlightlist = NULL;
386 int emax;
387 TRITREE *startleaf;
388 COLOR cval;
389 double growstep;
390 FVECT curcent;
391 double currad;
392 double curomega;
393 COLOR curintens;
394 double thisthresh;
395 int thisethresh;
396 int i;
397 /*
398 * General algorithm:
399 * 1) Start with brightest unclaimed pixel
400 * 2) Grow source toward brightest unclaimed perimeter until:
401 * a) Source exceeds maximum size, or
402 * b) Perimeter values all below threshold, or
403 * c) Source average drops below threshold
404 * 3) Loop until nothing over threshold
405 *
406 * Complexity added to absorb insignificant sources in larger ones.
407 */
408 if (thresh <= FTINY)
409 return;
410 for ( ; ; ) {
411 emax = ethresh; /* find brightest unclaimed */
412 startleaf = NULL;
413 for (i = 0; i < NTRUNKBR; i++) {
414 TRITREE *bigger = findemax(&samptree->kid[i], &emax);
415 if (bigger != NULL)
416 startleaf = bigger;
417 }
418 if (startleaf == NULL)
419 break;
420 /* claim it */
421 VCOPY(curcent, startleaf->sd);
422 curomega = tri_omegav(startleaf->gdv);
423 currad = sqrt(curomega/PI);
424 growstep = 3.*currad;
425 colr_color(curintens, startleaf->val);
426 thisthresh = .15*bright(curintens);
427 if (thisthresh < thresh) thisthresh = thresh;
428 thisethresh = (int)(log(thisthresh)/log(2.) + (COLXS+.5));
429 scalecolor(curintens, curomega);
430 copycolr(startleaf->val, blkclr);
431 do { /* grow source */
432 FVECT vsum;
433 double movedist;
434 vsum[0] = vsum[1] = vsum[2] = .0;
435 for (i = 0; i < NTRUNKBR; i++)
436 vector_sum(vsum, &samptree->kid[i],
437 curcent, cos(currad+growstep),
438 thisethresh);
439 if (normalize(vsum) == .0)
440 break;
441 movedist = acos(DOT(vsum,curcent));
442 if (movedist > growstep) {
443 VSUB(vsum, vsum, curcent);
444 movedist = growstep/VLEN(vsum);
445 VSUM(curcent, curcent, vsum, movedist);
446 normalize(curcent);
447 } else
448 VCOPY(curcent, vsum);
449 currad += growstep;
450 curomega = 2.*PI*(1. - cos(currad));
451 for (i = 0; i < NTRUNKBR; i++)
452 claimlight(curintens, &samptree->kid[i],
453 curcent, cos(currad));
454 } while (curomega < maxomega &&
455 bright(curintens)/curomega > thisthresh);
456 if (bright(curintens) < minintens) {
457 add2lost(&lostlightlist, curintens, curcent);
458 continue;
459 }
460 /* absorb lost contributions */
461 getlost(&lostlightlist, curintens, curcent, curomega);
462 ++nsrcs; /* output source */
463 scalecolor(curintens, 1./curomega);
464 printf("\nvoid illum IBLout\n");
465 printf("0\n0\n3 %f %f %f\n",
466 colval(curintens,RED),
467 colval(curintens,GRN),
468 colval(curintens,BLU));
469 printf("\nIBLout source IBLsrc%d\n", nsrcs);
470 printf("0\n0\n4 %f %f %f %f\n",
471 curcent[0], curcent[1], curcent[2],
472 2.*180./PI*currad);
473 }
474 }
475
476 int
477 main(int argc, char *argv[])
478 {
479 long nsamps = DEF_NSAMPS;
480 double maxang = DEF_MAXANG;
481 TRITREE *samptree;
482 double thresh = 0;
483 int i;
484
485 progname = argv[0];
486 for (i = 1; i < argc && argv[i][0] == '-'; i++)
487 switch (argv[i][1]) {
488 case 'd': /* number of samples */
489 if (i >= argc-1) goto userr;
490 nsamps = atol(argv[++i]);
491 break;
492 case 't': /* manual threshold */
493 if (i >= argc-1) goto userr;
494 thresh = atof(argv[++i]);
495 break;
496 case 'a': /* maximum source angle */
497 if (i >= argc-1) goto userr;
498 maxang = atof(argv[++i]);
499 if (maxang <= FTINY)
500 goto userr;
501 if (maxang > 180.)
502 maxang = 180.;
503 break;
504 default:
505 goto userr;
506 }
507 if (i < argc-1)
508 goto userr;
509 /* start our ray calculation */
510 directvis = 0;
511 ray_init(i == argc-1 ? argv[i] : (char *)NULL);
512 VCOPY(scene_cent, thescene.cuorg);
513 scene_cent[0] += 0.5*thescene.cusize;
514 scene_cent[1] += 0.5*thescene.cusize;
515 scene_cent[2] += 0.5*thescene.cusize;
516 scene_rad = 0.86603*thescene.cusize;
517 /* sample geodesic mesh */
518 samptree = geosample(nsamps);
519 /* get source threshold */
520 if (thresh <= FTINY)
521 thresh = get_threshold(samptree);
522 /* done with ray samples */
523 ray_done(1);
524 /* print header */
525 printf("# ");
526 printargs(argc, argv, stdout);
527 /* create & print sources */
528 mksources(samptree, thresh, maxang);
529 /* all done, no need to clean up */
530 return(0);
531 userr:
532 fprintf(stderr, "Usage: %s [-d nsamps][-t thresh][-a maxang] [octree]\n",
533 argv[0]);
534 exit(1);
535 }
536
537 void
538 eputs(char *s)
539 {
540 static int midline = 0;
541
542 if (!*s)
543 return;
544 if (!midline++) {
545 fputs(progname, stderr);
546 fputs(": ", stderr);
547 }
548 fputs(s, stderr);
549 if (s[strlen(s)-1] == '\n') {
550 fflush(stderr);
551 midline = 0;
552 }
553 }
554
555 void
556 wputs(char *s)
557 {
558 /* no warnings */
559 }