ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum2.c
Revision: 2.28
Committed: Thu Dec 13 07:03:37 2007 UTC (16 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R9
Changes since 2.27: +5 -3 lines
Log Message:
Superficially working version of mkillum with BSDF input

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mkillum2.c,v 2.27 2007/12/12 05:09:58 greg Exp $";
3 #endif
4 /*
5 * Routines to do the actual calculation for mkillum
6 */
7
8 #include <string.h>
9
10 #include "mkillum.h"
11 #include "face.h"
12 #include "cone.h"
13 #include "source.h"
14
15
16 COLORV * distarr = NULL; /* distribution array */
17 int distsiz = 0;
18
19
20 void
21 newdist( /* allocate & clear distribution array */
22 int siz
23 )
24 {
25 if (siz <= 0) {
26 if (distsiz > 0)
27 free((void *)distarr);
28 distarr = NULL;
29 distsiz = 0;
30 return;
31 }
32 if (distsiz < siz) {
33 if (distsiz > 0)
34 free((void *)distarr);
35 distarr = (COLORV *)malloc(sizeof(COLOR)*siz);
36 if (distarr == NULL)
37 error(SYSTEM, "out of memory in newdist");
38 distsiz = siz;
39 }
40 memset(distarr, '\0', sizeof(COLOR)*siz);
41 }
42
43
44 int
45 process_ray( /* process a ray result or report error */
46 RAY *r,
47 int rv
48 )
49 {
50 COLORV *colp;
51
52 if (rv == 0) /* no result ready */
53 return(0);
54 if (rv < 0)
55 error(USER, "ray tracing process died");
56 if (r->rno >= distsiz)
57 error(INTERNAL, "bad returned index in process_ray");
58 multcolor(r->rcol, r->rcoef); /* in case it's a source ray */
59 colp = &distarr[r->rno * 3];
60 addcolor(colp, r->rcol);
61 return(1);
62 }
63
64
65 void
66 raysamp( /* queue a ray sample */
67 int ndx,
68 FVECT org,
69 FVECT dir
70 )
71 {
72 RAY myRay;
73 int rv;
74
75 if ((ndx < 0) | (ndx >= distsiz))
76 error(INTERNAL, "bad index in raysamp");
77 VCOPY(myRay.rorg, org);
78 VCOPY(myRay.rdir, dir);
79 myRay.rmax = .0;
80 rayorigin(&myRay, PRIMARY, NULL, NULL);
81 myRay.rno = ndx;
82 /* queue ray, check result */
83 process_ray(&myRay, ray_pqueue(&myRay));
84 }
85
86
87 void
88 srcsamps( /* sample sources from this surface position */
89 struct illum_args *il,
90 FVECT org,
91 FVECT nrm,
92 MAT4 ixfm
93 )
94 {
95 int nalt, nazi;
96 SRCINDEX si;
97 RAY sr;
98 FVECT v;
99 double d;
100 int i, j;
101 /* get sampling density */
102 if (il->sampdens <= 0) {
103 nalt = nazi = 1;
104 } else {
105 i = PI * il->sampdens;
106 nalt = sqrt(i/PI) + .5;
107 nazi = PI*nalt + .5;
108 }
109 initsrcindex(&si); /* loop over (sub)sources */
110 for ( ; ; ) {
111 VCOPY(sr.rorg, org); /* pick side to shoot from */
112 if (il->sd != NULL) {
113 int sn = si.sn;
114 if (si.sp+1 >= si.np) ++sn;
115 if (sn >= nsources) break;
116 if (source[sn].sflags & SDISTANT)
117 d = DOT(source[sn].sloc, nrm);
118 else {
119 VSUB(v, source[sn].sloc, org);
120 d = DOT(v, nrm);
121 }
122 } else
123 d = 1.0; /* only transmission */
124 if (d < 0.0)
125 d = -1.0001*il->thick - 5.*FTINY;
126 else
127 d = 5.*FTINY;
128 for (i = 3; i--; )
129 sr.rorg[i] += d*nrm[i];
130 if (!srcray(&sr, NULL, &si))
131 break; /* end of sources */
132 /* index direction */
133 if (ixfm != NULL)
134 multv3(v, sr.rdir, ixfm);
135 else
136 VCOPY(v, sr.rdir);
137 if (il->sd != NULL) {
138 i = getBSDF_incndx(il->sd, v);
139 if (i < 0)
140 continue; /* must not be important */
141 sr.rno = i;
142 d = 1.0/getBSDF_incohm(il->sd, i);
143 } else {
144 if (v[2] >= -FTINY)
145 continue; /* only sample transmission */
146 d = 1.0 - v[2]*v[2];
147 i = d*nalt;
148 d = atan2(-v[1], -v[0])/(2.*PI);
149 if (d < 0.0) d += 1.0;
150 j = d*nazi + 0.5;
151 if (j >= nazi) j = 0;
152 sr.rno = i*nazi + j;
153 d = nalt*nazi/PI * -v[2];
154 }
155 d *= si.dom; /* solid angle correction */
156 scalecolor(sr.rcoef, d);
157 process_ray(&sr, ray_pqueue(&sr));
158 }
159 }
160
161
162 void
163 rayclean() /* finish all pending rays */
164 {
165 RAY myRay;
166
167 while (process_ray(&myRay, ray_presult(&myRay, 0)))
168 ;
169 }
170
171
172 static void
173 mkaxes( /* compute u and v to go with n */
174 FVECT u,
175 FVECT v,
176 FVECT n
177 )
178 {
179 register int i;
180
181 v[0] = v[1] = v[2] = 0.0;
182 for (i = 0; i < 3; i++)
183 if (n[i] < 0.6 && n[i] > -0.6)
184 break;
185 v[i] = 1.0;
186 fcross(u, v, n);
187 normalize(u);
188 fcross(v, n, u);
189 }
190
191
192 static void
193 rounddir( /* compute uniform spherical direction */
194 register FVECT dv,
195 double alt,
196 double azi
197 )
198 {
199 double d1, d2;
200
201 dv[2] = 1. - 2.*alt;
202 d1 = sqrt(1. - dv[2]*dv[2]);
203 d2 = 2.*PI * azi;
204 dv[0] = d1*cos(d2);
205 dv[1] = d1*sin(d2);
206 }
207
208
209 void
210 flatdir( /* compute uniform hemispherical direction */
211 FVECT dv,
212 double alt,
213 double azi
214 )
215 {
216 double d1, d2;
217
218 d1 = sqrt(alt);
219 d2 = 2.*PI * azi;
220 dv[0] = d1*cos(d2);
221 dv[1] = d1*sin(d2);
222 dv[2] = sqrt(1. - alt);
223 }
224
225
226 int
227 my_default( /* default illum action */
228 OBJREC *ob,
229 struct illum_args *il,
230 char *nm
231 )
232 {
233 sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
234 nm, ofun[ob->otype].funame, ob->oname);
235 error(WARNING, errmsg);
236 printobj(il->altmat, ob);
237 return(1);
238 }
239
240
241 int
242 my_face( /* make an illum face */
243 OBJREC *ob,
244 struct illum_args *il,
245 char *nm
246 )
247 {
248 int dim[2];
249 int n, nalt, nazi, alti;
250 double sp[2], r1, r2;
251 int h;
252 FVECT dn, org, dir;
253 FVECT u, v;
254 double ur[2], vr[2];
255 MAT4 xfm;
256 int nallow;
257 FACE *fa;
258 int i, j;
259 /* get/check arguments */
260 fa = getface(ob);
261 if (fa->area == 0.0) {
262 freeface(ob);
263 return(my_default(ob, il, nm));
264 }
265 /* set up sampling */
266 if (il->sd != NULL) {
267 if (!getBSDF_xfm(xfm, fa->norm, il->udir)) {
268 objerror(ob, WARNING, "illegal up direction");
269 freeface(ob);
270 return(my_default(ob, il, nm));
271 }
272 n = il->sd->ninc;
273 } else {
274 if (il->sampdens <= 0) {
275 nalt = nazi = 1; /* diffuse assumption */
276 } else {
277 n = PI * il->sampdens;
278 nalt = sqrt(n/PI) + .5;
279 nazi = PI*nalt + .5;
280 }
281 n = nazi*nalt;
282 }
283 newdist(n);
284 /* take first edge >= sqrt(area) */
285 for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
286 u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
287 u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
288 u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
289 if ((r1 = DOT(u,u)) >= fa->area-FTINY)
290 break;
291 }
292 if (i < fa->nv) { /* got one! -- let's align our axes */
293 r2 = 1.0/sqrt(r1);
294 u[0] *= r2; u[1] *= r2; u[2] *= r2;
295 fcross(v, fa->norm, u);
296 } else /* oh well, we'll just have to wing it */
297 mkaxes(u, v, fa->norm);
298 /* now, find limits in (u,v) coordinates */
299 ur[0] = vr[0] = FHUGE;
300 ur[1] = vr[1] = -FHUGE;
301 for (i = 0; i < fa->nv; i++) {
302 r1 = DOT(VERTEX(fa,i),u);
303 if (r1 < ur[0]) ur[0] = r1;
304 if (r1 > ur[1]) ur[1] = r1;
305 r2 = DOT(VERTEX(fa,i),v);
306 if (r2 < vr[0]) vr[0] = r2;
307 if (r2 > vr[1]) vr[1] = r2;
308 }
309 dim[0] = random();
310 /* sample polygon */
311 nallow = 5*n*il->nsamps;
312 for (dim[1] = 0; dim[1] < n; dim[1]++)
313 for (i = 0; i < il->nsamps; i++) {
314 /* randomize direction */
315 h = ilhash(dim, 2) + i;
316 if (il->sd != NULL) {
317 r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
318 } else {
319 multisamp(sp, 2, urand(h));
320 alti = dim[1]/nazi;
321 r1 = (alti + sp[0])/nalt;
322 r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
323 flatdir(dn, r1, r2);
324 for (j = 0; j < 3; j++)
325 dir[j] = -dn[0]*u[j] - dn[1]*v[j] -
326 dn[2]*fa->norm[j];
327 }
328 /* randomize location */
329 do {
330 multisamp(sp, 2, urand(h+4862+nallow));
331 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
332 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
333 for (j = 0; j < 3; j++)
334 org[j] = r1*u[j] + r2*v[j]
335 + fa->offset*fa->norm[j];
336 } while (!inface(org, fa) && nallow-- > 0);
337 if (nallow < 0) {
338 objerror(ob, WARNING, "bad aspect");
339 rayclean();
340 freeface(ob);
341 return(my_default(ob, il, nm));
342 }
343 if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY)
344 r1 = -1.0001*il->thick - 5.*FTINY;
345 else
346 r1 = 5.*FTINY;
347 for (j = 0; j < 3; j++)
348 org[j] += r1*fa->norm[j];
349 /* send sample */
350 raysamp(dim[1], org, dir);
351 }
352 /* add in direct component? */
353 if (!directvis && il->flags & IL_LIGHT) {
354 MAT4 ixfm;
355 if (il->sd == NULL) {
356 for (i = 3; i--; ) {
357 ixfm[i][0] = u[i];
358 ixfm[i][1] = v[i];
359 ixfm[i][2] = fa->norm[i];
360 ixfm[i][3] = 0.;
361 }
362 ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
363 ixfm[3][3] = 1.;
364 } else if (!invmat4(ixfm, xfm))
365 objerror(ob, INTERNAL, "cannot invert BSDF transform");
366 dim[0] = random();
367 nallow = 10*il->nsamps;
368 for (i = 0; i < il->nsamps; i++) {
369 /* randomize location */
370 h = dim[0] + samplendx++;
371 do {
372 multisamp(sp, 2, urand(h+nallow));
373 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
374 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
375 for (j = 0; j < 3; j++)
376 org[j] = r1*u[j] + r2*v[j]
377 + fa->offset*fa->norm[j];
378 } while (!inface(org, fa) && nallow-- > 0);
379 if (nallow < 0) {
380 objerror(ob, WARNING, "bad aspect");
381 rayclean();
382 freeface(ob);
383 return(my_default(ob, il, nm));
384 }
385 /* sample source rays */
386 srcsamps(il, org, fa->norm, ixfm);
387 }
388 }
389 /* wait for all rays to finish */
390 rayclean();
391 if (il->sd != NULL) { /* run distribution through BSDF */
392 nalt = sqrt(il->sd->nout/PI) + .5;
393 nazi = PI*nalt + .5;
394 redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm);
395 if (!il->sampdens)
396 il->sampdens = nalt*nazi/PI + .999;
397 }
398 /* write out the face and its distribution */
399 if (average(il, distarr, n)) {
400 if (il->sampdens > 0)
401 flatout(il, distarr, nalt, nazi, u, v, fa->norm);
402 illumout(il, ob);
403 } else
404 printobj(il->altmat, ob);
405 /* clean up */
406 freeface(ob);
407 return(0);
408 }
409
410
411 int
412 my_sphere( /* make an illum sphere */
413 register OBJREC *ob,
414 struct illum_args *il,
415 char *nm
416 )
417 {
418 int dim[3];
419 int n, nalt, nazi;
420 double sp[4], r1, r2, r3;
421 FVECT org, dir;
422 FVECT u, v;
423 register int i, j;
424 /* check arguments */
425 if (ob->oargs.nfargs != 4)
426 objerror(ob, USER, "bad # of arguments");
427 /* set up sampling */
428 if (il->sampdens <= 0)
429 nalt = nazi = 1;
430 else {
431 n = 4.*PI * il->sampdens;
432 nalt = sqrt(2./PI*n) + .5;
433 nazi = PI/2.*nalt + .5;
434 }
435 if (il->sd != NULL)
436 objerror(ob, WARNING, "BSDF ignored");
437 n = nalt*nazi;
438 newdist(n);
439 dim[0] = random();
440 /* sample sphere */
441 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
442 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
443 for (i = 0; i < il->nsamps; i++) {
444 /* next sample point */
445 multisamp(sp, 4, urand(ilhash(dim,3)+i));
446 /* random direction */
447 r1 = (dim[1] + sp[0])/nalt;
448 r2 = (dim[2] + sp[1] - .5)/nazi;
449 rounddir(dir, r1, r2);
450 /* random location */
451 mkaxes(u, v, dir); /* yuck! */
452 r3 = sqrt(sp[2]);
453 r2 = 2.*PI*sp[3];
454 r1 = r3*ob->oargs.farg[3]*cos(r2);
455 r2 = r3*ob->oargs.farg[3]*sin(r2);
456 r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
457 for (j = 0; j < 3; j++) {
458 org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
459 r3*dir[j];
460 dir[j] = -dir[j];
461 }
462 /* send sample */
463 raysamp(dim[1]*nazi+dim[2], org, dir);
464 }
465 /* wait for all rays to finish */
466 rayclean();
467 /* write out the sphere and its distribution */
468 if (average(il, distarr, n)) {
469 if (il->sampdens > 0)
470 roundout(il, distarr, nalt, nazi);
471 else
472 objerror(ob, WARNING, "diffuse distribution");
473 illumout(il, ob);
474 } else
475 printobj(il->altmat, ob);
476 /* clean up */
477 return(1);
478 }
479
480
481 int
482 my_ring( /* make an illum ring */
483 OBJREC *ob,
484 struct illum_args *il,
485 char *nm
486 )
487 {
488 int dim[2];
489 int n, nalt, nazi, alti;
490 double sp[2], r1, r2, r3;
491 int h;
492 FVECT dn, org, dir;
493 FVECT u, v;
494 MAT4 xfm;
495 CONE *co;
496 int i, j;
497 /* get/check arguments */
498 co = getcone(ob, 0);
499 /* set up sampling */
500 if (il->sd != NULL) {
501 if (!getBSDF_xfm(xfm, co->ad, il->udir)) {
502 objerror(ob, WARNING, "illegal up direction");
503 freecone(ob);
504 return(my_default(ob, il, nm));
505 }
506 n = il->sd->ninc;
507 } else {
508 if (il->sampdens <= 0) {
509 nalt = nazi = 1; /* diffuse assumption */
510 } else {
511 n = PI * il->sampdens;
512 nalt = sqrt(n/PI) + .5;
513 nazi = PI*nalt + .5;
514 }
515 n = nazi*nalt;
516 }
517 newdist(n);
518 mkaxes(u, v, co->ad);
519 dim[0] = random();
520 /* sample disk */
521 for (dim[1] = 0; dim[1] < n; dim[1]++)
522 for (i = 0; i < il->nsamps; i++) {
523 /* next sample point */
524 h = ilhash(dim,2) + i;
525 /* randomize direction */
526 if (il->sd != NULL) {
527 r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
528 } else {
529 multisamp(sp, 2, urand(h));
530 alti = dim[1]/nazi;
531 r1 = (alti + sp[0])/nalt;
532 r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
533 flatdir(dn, r1, r2);
534 for (j = 0; j < 3; j++)
535 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
536 }
537 /* randomize location */
538 multisamp(sp, 2, urand(h+8371));
539 r3 = sqrt(CO_R0(co)*CO_R0(co) +
540 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
541 r2 = 2.*PI*sp[1];
542 r1 = r3*cos(r2);
543 r2 = r3*sin(r2);
544 if (il->sd != NULL && DOT(dir, co->ad) < -FTINY)
545 r3 = -1.0001*il->thick - 5.*FTINY;
546 else
547 r3 = 5.*FTINY;
548 for (j = 0; j < 3; j++)
549 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
550 r3*co->ad[j];
551 /* send sample */
552 raysamp(dim[1], org, dir);
553 }
554 /* add in direct component? */
555 if (!directvis && il->flags & IL_LIGHT) {
556 MAT4 ixfm;
557 if (il->sd == NULL) {
558 for (i = 3; i--; ) {
559 ixfm[i][0] = u[i];
560 ixfm[i][1] = v[i];
561 ixfm[i][2] = co->ad[i];
562 ixfm[i][3] = 0.;
563 }
564 ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
565 ixfm[3][3] = 1.;
566 } else if (!invmat4(ixfm, xfm))
567 objerror(ob, INTERNAL, "cannot invert BSDF transform");
568 dim[0] = random();
569 for (i = 0; i < il->nsamps; i++) {
570 /* randomize location */
571 h = dim[0] + samplendx++;
572 multisamp(sp, 2, urand(h));
573 r3 = sqrt(CO_R0(co)*CO_R0(co) +
574 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
575 r2 = 2.*PI*sp[1];
576 r1 = r3*cos(r2);
577 r2 = r3*sin(r2);
578 for (j = 0; j < 3; j++)
579 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j];
580 /* sample source rays */
581 srcsamps(il, org, co->ad, ixfm);
582 }
583 }
584 /* wait for all rays to finish */
585 rayclean();
586 if (il->sd != NULL) { /* run distribution through BSDF */
587 nalt = sqrt(il->sd->nout/PI) + .5;
588 nazi = PI*nalt + .5;
589 redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
590 if (!il->sampdens)
591 il->sampdens = nalt*nazi/PI + .999;
592 }
593 /* write out the ring and its distribution */
594 if (average(il, distarr, n)) {
595 if (il->sampdens > 0)
596 flatout(il, distarr, nalt, nazi, u, v, co->ad);
597 illumout(il, ob);
598 } else
599 printobj(il->altmat, ob);
600 /* clean up */
601 freecone(ob);
602 return(1);
603 }