ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mksource.c
Revision: 2.3
Committed: Mon Sep 19 12:32:12 2005 UTC (18 years, 7 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.2: +3 -2 lines
Log Message:
Added mksource to SCons build, fixed some compile warnings.

File Contents

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