ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum2.c
Revision: 2.26
Committed: Sat Dec 8 01:43:09 2007 UTC (16 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.25: +2 -3 lines
Log Message:
Additional changes towards BSDF implementation

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mkillum2.c,v 2.25 2007/12/06 19:40:43 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;
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 }
396 /* write out the face and its distribution */
397 if (average(il, distarr, n)) {
398 if (il->sampdens > 0)
399 flatout(il, distarr, nalt, nazi, u, v, fa->norm);
400 illumout(il, ob);
401 } else
402 printobj(il->altmat, ob);
403 /* clean up */
404 freeface(ob);
405 return(0);
406 }
407
408
409 int
410 my_sphere( /* make an illum sphere */
411 register OBJREC *ob,
412 struct illum_args *il,
413 char *nm
414 )
415 {
416 int dim[3];
417 int n, nalt, nazi;
418 double sp[4], r1, r2, r3;
419 FVECT org, dir;
420 FVECT u, v;
421 register int i, j;
422 /* check arguments */
423 if (ob->oargs.nfargs != 4)
424 objerror(ob, USER, "bad # of arguments");
425 /* set up sampling */
426 if (il->sampdens <= 0)
427 nalt = nazi = 1;
428 else {
429 n = 4.*PI * il->sampdens;
430 nalt = sqrt(2./PI*n) + .5;
431 nazi = PI/2.*nalt + .5;
432 }
433 if (il->sd != NULL)
434 objerror(ob, WARNING, "BSDF ignored");
435 n = nalt*nazi;
436 newdist(n);
437 dim[0] = random();
438 /* sample sphere */
439 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
440 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
441 for (i = 0; i < il->nsamps; i++) {
442 /* next sample point */
443 multisamp(sp, 4, urand(ilhash(dim,3)+i));
444 /* random direction */
445 r1 = (dim[1] + sp[0])/nalt;
446 r2 = (dim[2] + sp[1] - .5)/nazi;
447 rounddir(dir, r1, r2);
448 /* random location */
449 mkaxes(u, v, dir); /* yuck! */
450 r3 = sqrt(sp[2]);
451 r2 = 2.*PI*sp[3];
452 r1 = r3*ob->oargs.farg[3]*cos(r2);
453 r2 = r3*ob->oargs.farg[3]*sin(r2);
454 r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
455 for (j = 0; j < 3; j++) {
456 org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
457 r3*dir[j];
458 dir[j] = -dir[j];
459 }
460 /* send sample */
461 raysamp(dim[1]*nazi+dim[2], org, dir);
462 }
463 /* wait for all rays to finish */
464 rayclean();
465 /* write out the sphere and its distribution */
466 if (average(il, distarr, n)) {
467 if (il->sampdens > 0)
468 roundout(il, distarr, nalt, nazi);
469 else
470 objerror(ob, WARNING, "diffuse distribution");
471 illumout(il, ob);
472 } else
473 printobj(il->altmat, ob);
474 /* clean up */
475 return(1);
476 }
477
478
479 int
480 my_ring( /* make an illum ring */
481 OBJREC *ob,
482 struct illum_args *il,
483 char *nm
484 )
485 {
486 int dim[2];
487 int n, nalt, nazi, alti;
488 double sp[2], r1, r2, r3;
489 int h;
490 FVECT dn, org, dir;
491 FVECT u, v;
492 MAT4 xfm;
493 CONE *co;
494 int i, j;
495 /* get/check arguments */
496 co = getcone(ob, 0);
497 /* set up sampling */
498 if (il->sd != NULL) {
499 if (!getBSDF_xfm(xfm, co->ad, il->udir)) {
500 objerror(ob, WARNING, "illegal up direction");
501 freecone(ob);
502 return(my_default(ob, il, nm));
503 }
504 n = il->sd->ninc;
505 } else {
506 if (il->sampdens <= 0) {
507 nalt = nazi = 1; /* diffuse assumption */
508 } else {
509 n = PI * il->sampdens;
510 nalt = sqrt(n/PI) + .5;
511 nazi = PI*nalt + .5;
512 }
513 n = nazi*nalt;
514 }
515 newdist(n);
516 mkaxes(u, v, co->ad);
517 dim[0] = random();
518 /* sample disk */
519 for (dim[1] = 0; dim[1] < n; dim[1]++)
520 for (i = 0; i < il->nsamps; i++) {
521 /* next sample point */
522 h = ilhash(dim,2) + i;
523 /* randomize direction */
524 if (il->sd != NULL) {
525 r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
526 } else {
527 multisamp(sp, 2, urand(h));
528 alti = dim[1]/nazi;
529 r1 = (alti + sp[0])/nalt;
530 r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
531 flatdir(dn, r1, r2);
532 for (j = 0; j < 3; j++)
533 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
534 }
535 /* randomize location */
536 multisamp(sp, 2, urand(h+8371));
537 r3 = sqrt(CO_R0(co)*CO_R0(co) +
538 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
539 r2 = 2.*PI*sp[1];
540 r1 = r3*cos(r2);
541 r2 = r3*sin(r2);
542 if (il->sd != NULL && DOT(dir, co->ad) < -FTINY)
543 r3 = -1.0001*il->thick - 5.*FTINY;
544 else
545 r3 = 5.*FTINY;
546 for (j = 0; j < 3; j++)
547 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
548 r3*co->ad[j];
549 /* send sample */
550 raysamp(dim[1], org, dir);
551 }
552 /* add in direct component? */
553 if (!directvis && il->flags & IL_LIGHT) {
554 MAT4 ixfm;
555 if (il->sd == NULL) {
556 for (i = 3; i--; ) {
557 ixfm[i][0] = u[i];
558 ixfm[i][1] = v[i];
559 ixfm[i][2] = co->ad[i];
560 ixfm[i][3] = 0.;
561 }
562 ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
563 ixfm[3][3] = 1.;
564 } else if (!invmat4(ixfm, xfm))
565 objerror(ob, INTERNAL, "cannot invert BSDF transform");
566 dim[0] = random();
567 for (i = 0; i < il->nsamps; i++) {
568 /* randomize location */
569 h = dim[0] + samplendx++;
570 multisamp(sp, 2, urand(h));
571 r3 = sqrt(CO_R0(co)*CO_R0(co) +
572 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
573 r2 = 2.*PI*sp[1];
574 r1 = r3*cos(r2);
575 r2 = r3*sin(r2);
576 for (j = 0; j < 3; j++)
577 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j];
578 /* sample source rays */
579 srcsamps(il, org, co->ad, ixfm);
580 }
581 }
582 /* wait for all rays to finish */
583 rayclean();
584 if (il->sd != NULL) { /* run distribution through BSDF */
585 nalt = sqrt(il->sd->nout/PI) + .5;
586 nazi = PI*nalt + .5;
587 redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
588 }
589 /* write out the ring and its distribution */
590 if (average(il, distarr, n)) {
591 if (il->sampdens > 0)
592 flatout(il, distarr, nalt, nazi, u, v, co->ad);
593 illumout(il, ob);
594 } else
595 printobj(il->altmat, ob);
596 /* clean up */
597 freecone(ob);
598 return(1);
599 }