ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mksource.c
Revision: 2.5
Committed: Fri Aug 24 04:41:36 2007 UTC (16 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0, rad3R9
Changes since 2.4: +89 -71 lines
Log Message:
Made ray direction storage more efficient (1/4 the original memory)

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mksource.c,v 2.4 2005/09/23 19:22:37 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 setcolr(leaf->val, colval(myray.rcol,RED),
120 colval(myray.rcol,GRN),
121 colval(myray.rcol,BLU));
122 }
123
124 /* Initialize a branch node contained in the given spherical triangle */
125 void
126 subdivide(TRITREE *branch, const int32 dv[3])
127 {
128 FVECT dvv[3], sdv[3];
129 int32 sd[3];
130 int i;
131
132 for (i = 0; i < 3; i++) { /* copy spherical triangle */
133 branch->gdv[i] = dv[i];
134 decodedir(dvv[i], dv[i]);
135 }
136 for (i = 0; i < 3; i++) { /* create new vertices */
137 int j = (i+1)%3;
138 VADD(sdv[i], dvv[i], dvv[j]);
139 normalize(sdv[i]);
140 sd[i] = encodedir(sdv[i]);
141 }
142 /* allocate leaves */
143 branch->kid = (TRITREE *)calloc(4, sizeof(TRITREE));
144 if (branch->kid == NULL)
145 error(SYSTEM, "out of memory in subdivide()");
146 /* assign subtriangle directions */
147 branch->kid[0].gdv[0] = dv[0];
148 branch->kid[0].gdv[1] = sd[0];
149 branch->kid[0].gdv[2] = sd[2];
150 branch->kid[1].gdv[0] = sd[0];
151 branch->kid[1].gdv[1] = dv[1];
152 branch->kid[1].gdv[2] = sd[1];
153 branch->kid[2].gdv[0] = sd[1];
154 branch->kid[2].gdv[1] = dv[2];
155 branch->kid[2].gdv[2] = sd[2];
156 branch->kid[3].gdv[0] = sd[0];
157 branch->kid[3].gdv[1] = sd[1];
158 branch->kid[3].gdv[2] = sd[2];
159 }
160
161 /* Recursively subdivide the given node to the specified quadtree depth */
162 void
163 branchsample(TRITREE *node, int depth)
164 {
165 int i;
166
167 if (depth <= 0)
168 return;
169 if (isleaf(node)) { /* subdivide leaf node */
170 TRITREE branch, *moved_leaf;
171 FVECT sdir;
172 subdivide(&branch, node->gdv);
173 decodedir(sdir, node->sd);
174 moved_leaf = findleaf(&branch, sdir);
175 if (moved_leaf != NULL) { /* bequeath old sample */
176 moved_leaf->sd = node->sd;
177 copycolr(moved_leaf->val, node->val);
178 }
179 for (i = 0; i < 4; i++) /* compute new samples */
180 if (&branch.kid[i] != moved_leaf)
181 leafsample(&branch.kid[i]);
182 *node = branch; /* replace leaf with branch */
183 }
184 for (i = 0; i < 4; i++) /* subdivide children */
185 branchsample(&node->kid[i], depth-1);
186 }
187
188 /* Sample sphere using triangular geodesic mesh */
189 TRITREE *
190 geosample(int nsamps)
191 {
192 int depth;
193 TRITREE *tree;
194 FVECT trunk[NTRUNKVERT];
195 int i, j;
196 /* figure out depth */
197 if ((nsamps -= 4) < 0)
198 error(USER, "minimum number of samples is 4");
199 nsamps = nsamps*3/NTRUNKBR; /* round up */
200 for (depth = 0; nsamps > 1; depth++)
201 nsamps >>= 2;
202 /* make base tetrahedron */
203 tree = (TRITREE *)malloc(sizeof(TRITREE));
204 if (tree == NULL) goto memerr;
205 trunk[0][0] = trunk[0][1] = 0; trunk[0][2] = 1;
206 trunk[1][0] = 0;
207 trunk[1][2] = cos(2.*asin(sqrt(2./3.)));
208 trunk[1][1] = sqrt(1. - trunk[1][2]*trunk[1][2]);
209 spinvector(trunk[2], trunk[1], trunk[0], 2.*PI/3.);
210 spinvector(trunk[3], trunk[1], trunk[0], 4.*PI/3.);
211 tree->gdv[0] = tree->gdv[1] = tree->gdv[2] = encodedir(trunk[0]);
212 tree->kid = (TRITREE *)calloc(NTRUNKBR, sizeof(TRITREE));
213 if (tree->kid == NULL) goto memerr;
214 /* grow our tree from trunk */
215 for (i = 0; i < NTRUNKBR; i++) {
216 for (j = 0; j < 3; j++) /* XXX works for tetra only */
217 tree->kid[i].gdv[j] = encodedir(trunk[(i+j)%NTRUNKVERT]);
218 leafsample(&tree->kid[i]);
219 branchsample(&tree->kid[i], depth);
220 }
221 return(tree);
222 memerr:
223 error(SYSTEM, "out of memory in geosample()");
224 return NULL; /* dummy return */
225 }
226
227 /* Compute leaf exponent histogram */
228 void
229 get_ehisto(const TRITREE *node, long exphisto[256])
230 {
231 int i;
232
233 if (isleaf(node)) {
234 ++exphisto[node->val[EXP]];
235 return;
236 }
237 for (i = 0; i < 4; i++)
238 get_ehisto(&node->kid[i], exphisto);
239 }
240
241 /* Get reasonable source threshold */
242 double
243 get_threshold(const TRITREE *tree)
244 {
245 long exphisto[256];
246 long samptotal;
247 int i;
248 /* compute sample histogram */
249 memset((void *)exphisto, 0, sizeof(exphisto));
250 for (i = 0; i < NTRUNKBR; i++)
251 get_ehisto(&tree->kid[i], exphisto);
252 /* use 98th percentile */
253 for (i = 0; i < 256; i++)
254 samptotal += exphisto[i];
255 samptotal /= 50;
256 for (i = 256; (--i > 0) & (samptotal > 0); )
257 samptotal -= exphisto[i];
258 return(ldexp(.75, i-COLXS));
259 }
260
261 /* Find leaf containing the maximum exponent */
262 TRITREE *
263 findemax(TRITREE *node, int *expp)
264 {
265 if (!isleaf(node)) {
266 TRITREE *maxleaf;
267 TRITREE *rleaf;
268 maxleaf = findemax(&node->kid[0], expp);
269 rleaf = findemax(&node->kid[1], expp);
270 if (rleaf != NULL) maxleaf = rleaf;
271 rleaf = findemax(&node->kid[2], expp);
272 if (rleaf != NULL) maxleaf = rleaf;
273 rleaf = findemax(&node->kid[3], expp);
274 if (rleaf != NULL) maxleaf = rleaf;
275 return(maxleaf);
276 }
277 if (node->val[EXP] <= *expp)
278 return(NULL);
279 *expp = node->val[EXP];
280 return(node);
281 }
282
283 /* Compute solid angle of spherical triangle (approx.) */
284 double
285 tri_omegav(const int32 vd[3])
286 {
287 FVECT v[3], e1, e2, vcross;
288
289 decodedir(v[0], vd[0]);
290 decodedir(v[1], vd[1]);
291 decodedir(v[2], vd[2]);
292 VSUB(e1, v[1], v[0]);
293 VSUB(e2, v[2], v[1]);
294 fcross(vcross, e1, e2);
295 return(.5*VLEN(vcross));
296 }
297
298 /* Sum intensity times direction for above-threshold perimiter within radius */
299 void
300 vector_sum(FVECT vsum, TRITREE *node,
301 FVECT cent, double maxr2, int ethresh)
302 {
303 if (isleaf(node)) {
304 double intens;
305 FVECT sdir;
306 if (node->val[EXP] < ethresh)
307 return; /* below threshold */
308 if (fdir2diff(node->sd,cent) > maxr2)
309 return; /* too far away */
310 intens = colrval(node->val,GRN) * tri_omegav(node->gdv);
311 decodedir(sdir, node->sd);
312 VSUM(vsum, vsum, sdir, intens);
313 return;
314 }
315 if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
316 fdir2diff(node->gdv[0],cent) < maxr2 &&
317 fdir2diff(node->gdv[1],cent) < maxr2 &&
318 fdir2diff(node->gdv[2],cent) < maxr2)
319 return; /* containing node */
320 vector_sum(vsum, &node->kid[0], cent, maxr2, ethresh);
321 vector_sum(vsum, &node->kid[1], cent, maxr2, ethresh);
322 vector_sum(vsum, &node->kid[2], cent, maxr2, ethresh);
323 vector_sum(vsum, &node->kid[3], cent, maxr2, ethresh);
324 }
325
326 /* Claim source contributions within the given solid angle */
327 void
328 claimlight(COLOR intens, TRITREE *node, FVECT cent, double maxr2)
329 {
330 int remaining;
331 int i;
332 if (isleaf(node)) { /* claim contribution */
333 COLOR contrib;
334 if (node->val[EXP] <= 0)
335 return; /* already claimed */
336 if (fdir2diff(node->sd,cent) > maxr2)
337 return; /* too far away */
338 colr_color(contrib, node->val);
339 scalecolor(contrib, tri_omegav(node->gdv));
340 addcolor(intens, contrib);
341 copycolr(node->val, blkclr);
342 return;
343 }
344 if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
345 fdir2diff(node->gdv[0],cent) < maxr2 &&
346 fdir2diff(node->gdv[1],cent) < maxr2 &&
347 fdir2diff(node->gdv[2],cent) < maxr2)
348 return; /* previously claimed node */
349 remaining = 0; /* recurse on children */
350 for (i = 0; i < 4; i++) {
351 claimlight(intens, &node->kid[i], cent, maxr2);
352 if (!isleaf(&node->kid[i]) || node->kid[i].val[EXP] != 0)
353 ++remaining;
354 }
355 if (remaining)
356 return;
357 /* consolidate empties */
358 free((void *)node->kid); node->kid = NULL;
359 copycolr(node->val, blkclr);
360 node->sd = node->gdv[0]; /* doesn't really matter */
361 }
362
363 /* Add lost light contribution to the given list */
364 void
365 add2lost(LOSTLIGHT **llp, COLOR intens, FVECT cent)
366 {
367 LOSTLIGHT *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
368
369 if (newll == NULL)
370 return;
371 copycolor(newll->intens, intens);
372 newll->sd = encodedir(cent);
373 newll->next = *llp;
374 *llp = newll;
375 }
376
377 /* Check lost light list for contributions */
378 void
379 getlost(LOSTLIGHT **llp, COLOR intens, FVECT cent, double omega)
380 {
381 const double maxr2 = omega/PI;
382 LOSTLIGHT lhead, *lastp, *thisp;
383
384 lhead.next = *llp;
385 lastp = &lhead;
386 while ((thisp = lastp->next) != NULL)
387 if (fdir2diff(thisp->sd,cent) <= maxr2) {
388 LOSTLIGHT *mynext = thisp->next;
389 addcolor(intens, thisp->intens);
390 free((void *)thisp);
391 lastp->next = mynext;
392 } else
393 lastp = thisp;
394 *llp = lhead.next;
395 }
396
397 /* Create & print distant sources */
398 void
399 mksources(TRITREE *samptree, double thresh, double maxang)
400 {
401 const int ethresh = (int)(log(thresh)/log(2.) + (COLXS+.5));
402 const double maxomega = 2.*PI*(1. - cos(PI/180./2.*maxang));
403 const double minintens = .05*thresh*maxomega;
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 for ( ; ; ) {
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 }
493 }
494
495 int
496 main(int argc, char *argv[])
497 {
498 long nsamps = DEF_NSAMPS;
499 double maxang = DEF_MAXANG;
500 TRITREE *samptree;
501 double thresh = 0;
502 int i;
503
504 progname = argv[0];
505 for (i = 1; i < argc && argv[i][0] == '-'; i++)
506 switch (argv[i][1]) {
507 case 'd': /* number of samples */
508 if (i >= argc-1) goto userr;
509 nsamps = atol(argv[++i]);
510 break;
511 case 't': /* manual threshold */
512 if (i >= argc-1) goto userr;
513 thresh = atof(argv[++i]);
514 break;
515 case 'a': /* maximum source angle */
516 if (i >= argc-1) goto userr;
517 maxang = atof(argv[++i]);
518 if (maxang <= FTINY)
519 goto userr;
520 if (maxang > 180.)
521 maxang = 180.;
522 break;
523 default:
524 goto userr;
525 }
526 if (i < argc-1)
527 goto userr;
528 /* start our ray calculation */
529 directvis = 0;
530 ray_init(i == argc-1 ? argv[i] : (char *)NULL);
531 VCOPY(scene_cent, thescene.cuorg);
532 scene_cent[0] += 0.5*thescene.cusize;
533 scene_cent[1] += 0.5*thescene.cusize;
534 scene_cent[2] += 0.5*thescene.cusize;
535 scene_rad = 0.86603*thescene.cusize;
536 /* sample geodesic mesh */
537 samptree = geosample(nsamps);
538 /* get source threshold */
539 if (thresh <= FTINY)
540 thresh = get_threshold(samptree);
541 /* done with ray samples */
542 ray_done(1);
543 /* print header */
544 printf("# ");
545 printargs(argc, argv, stdout);
546 /* create & print sources */
547 mksources(samptree, thresh, maxang);
548 /* all done, no need to clean up */
549 return(0);
550 userr:
551 fprintf(stderr, "Usage: %s [-d nsamps][-t thresh][-a maxang] [octree]\n",
552 argv[0]);
553 exit(1);
554 }
555
556 void
557 eputs(char *s)
558 {
559 static int midline = 0;
560
561 if (!*s)
562 return;
563 if (!midline++) {
564 fputs(progname, stderr);
565 fputs(": ", stderr);
566 }
567 fputs(s, stderr);
568 if (s[strlen(s)-1] == '\n') {
569 fflush(stderr);
570 midline = 0;
571 }
572 }
573
574 void
575 wputs(char *s)
576 {
577 /* no warnings */
578 }