ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum2.c
Revision: 2.34
Committed: Wed Sep 9 15:32:20 2009 UTC (14 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0
Changes since 2.33: +2 -2 lines
Log Message:
Fixed bug introduced in last bug-fix.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mkillum2.c,v 2.33 2009/09/08 23:05:47 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 #ifndef NBSDFSAMPS
16 #define NBSDFSAMPS 256 /* BSDF resampling count */
17 #endif
18
19 COLORV * distarr = NULL; /* distribution array */
20 int distsiz = 0;
21 COLORV * direct_discount = NULL; /* amount to take off direct */
22
23
24 void
25 newdist( /* allocate & clear distribution array */
26 int siz
27 )
28 {
29 if (siz <= 0) {
30 if (distsiz > 0)
31 free((void *)distarr);
32 distarr = NULL;
33 distsiz = 0;
34 return;
35 }
36 if (distsiz < siz) {
37 if (distsiz > 0)
38 free((void *)distarr);
39 distarr = (COLORV *)malloc(sizeof(COLOR)*siz);
40 if (distarr == NULL)
41 error(SYSTEM, "out of memory in newdist");
42 distsiz = siz;
43 }
44 memset(distarr, '\0', sizeof(COLOR)*siz);
45 }
46
47
48 static void
49 new_discount() /* allocate space for direct contrib. record */
50 {
51 if (distsiz <= 0)
52 return;
53 direct_discount = (COLORV *)calloc(distsiz, sizeof(COLOR));
54 if (direct_discount == NULL)
55 error(SYSTEM, "out of memory in new_discount");
56 }
57
58
59 static void
60 done_discount() /* clear off direct contrib. record */
61 {
62 if (direct_discount == NULL)
63 return;
64 free((void *)direct_discount);
65 direct_discount = NULL;
66 }
67
68
69 int
70 process_ray( /* process a ray result or report error */
71 RAY *r,
72 int rv
73 )
74 {
75 COLORV *colp;
76
77 if (rv == 0) /* no result ready */
78 return(0);
79 if (rv < 0)
80 error(USER, "ray tracing process died");
81 if (r->rno >= distsiz)
82 error(INTERNAL, "bad returned index in process_ray");
83 multcolor(r->rcol, r->rcoef); /* in case it's a source ray */
84 colp = &distarr[r->rno * 3];
85 addcolor(colp, r->rcol);
86 if (r->rsrc >= 0 && /* remember source contrib. */
87 direct_discount != NULL) {
88 colp = &direct_discount[r->rno * 3];
89 addcolor(colp, r->rcol);
90 }
91 return(1);
92 }
93
94
95 void
96 raysamp( /* queue a ray sample */
97 int ndx,
98 FVECT org,
99 FVECT dir
100 )
101 {
102 RAY myRay;
103 int rv;
104
105 if ((ndx < 0) | (ndx >= distsiz))
106 error(INTERNAL, "bad index in raysamp");
107 VCOPY(myRay.rorg, org);
108 VCOPY(myRay.rdir, dir);
109 myRay.rmax = .0;
110 rayorigin(&myRay, PRIMARY, NULL, NULL);
111 myRay.rno = ndx;
112 /* queue ray, check result */
113 process_ray(&myRay, ray_pqueue(&myRay));
114 }
115
116
117 void
118 srcsamps( /* sample sources from this surface position */
119 struct illum_args *il,
120 FVECT org,
121 FVECT nrm,
122 MAT4 ixfm
123 )
124 {
125 int nalt, nazi;
126 SRCINDEX si;
127 RAY sr;
128 FVECT v;
129 double d;
130 int i, j;
131 /* get sampling density */
132 if (il->sampdens <= 0) {
133 nalt = nazi = 1;
134 } else {
135 i = PI * il->sampdens;
136 nalt = sqrt(i/PI) + .5;
137 nazi = PI*nalt + .5;
138 }
139 initsrcindex(&si); /* loop over (sub)sources */
140 for ( ; ; ) {
141 VCOPY(sr.rorg, org); /* pick side to shoot from */
142 if (il->sd != NULL) {
143 int sn = si.sn;
144 if (si.sp+1 >= si.np) ++sn;
145 if (sn >= nsources) break;
146 if (source[sn].sflags & SDISTANT)
147 d = DOT(source[sn].sloc, nrm);
148 else {
149 VSUB(v, source[sn].sloc, org);
150 d = DOT(v, nrm);
151 }
152 } else
153 d = 1.0; /* only transmission */
154 if (d < 0.0)
155 d = -1.0001*il->thick - 5.*FTINY;
156 else
157 d = 5.*FTINY;
158 for (i = 3; i--; )
159 sr.rorg[i] += d*nrm[i];
160 samplendx++; /* increment sample counter */
161 if (!srcray(&sr, NULL, &si))
162 break; /* end of sources */
163 /* index direction */
164 if (ixfm != NULL)
165 multv3(v, sr.rdir, ixfm);
166 else
167 VCOPY(v, sr.rdir);
168 if (il->sd != NULL) {
169 i = getBSDF_incndx(il->sd, v);
170 if (i < 0)
171 continue; /* must not be important */
172 sr.rno = i;
173 d = 1.0/getBSDF_incohm(il->sd, i);
174 } else {
175 if (v[2] >= -FTINY)
176 continue; /* only sample transmission */
177 v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2];
178 sr.rno = flatindex(v, nalt, nazi);
179 d = nalt*nazi*(1./PI) * v[2];
180 }
181 d *= si.dom; /* solid angle correction */
182 scalecolor(sr.rcoef, d);
183 process_ray(&sr, ray_pqueue(&sr));
184 }
185 }
186
187
188 void
189 rayclean() /* finish all pending rays */
190 {
191 RAY myRay;
192
193 while (process_ray(&myRay, ray_presult(&myRay, 0)))
194 ;
195 }
196
197
198 static void
199 mkaxes( /* compute u and v to go with n */
200 FVECT u,
201 FVECT v,
202 FVECT n
203 )
204 {
205 register int i;
206
207 v[0] = v[1] = v[2] = 0.0;
208 for (i = 0; i < 3; i++)
209 if (n[i] < 0.6 && n[i] > -0.6)
210 break;
211 v[i] = 1.0;
212 fcross(u, v, n);
213 normalize(u);
214 fcross(v, n, u);
215 }
216
217
218 static void
219 rounddir( /* compute uniform spherical direction */
220 register FVECT dv,
221 double alt,
222 double azi
223 )
224 {
225 double d1, d2;
226
227 dv[2] = 1. - 2.*alt;
228 d1 = sqrt(1. - dv[2]*dv[2]);
229 d2 = 2.*PI * azi;
230 dv[0] = d1*cos(d2);
231 dv[1] = d1*sin(d2);
232 }
233
234
235 void
236 flatdir( /* compute uniform hemispherical direction */
237 FVECT dv,
238 double alt,
239 double azi
240 )
241 {
242 double d1, d2;
243
244 d1 = sqrt(alt);
245 d2 = 2.*PI * azi;
246 dv[0] = d1*cos(d2);
247 dv[1] = d1*sin(d2);
248 dv[2] = sqrt(1. - alt);
249 }
250
251 int
252 flatindex( /* compute index for hemispherical direction */
253 FVECT dv,
254 int nalt,
255 int nazi
256 )
257 {
258 double d;
259 int i, j;
260
261 d = 1.0 - dv[2]*dv[2];
262 i = d*nalt;
263 d = atan2(dv[1], dv[0]) * (0.5/PI);
264 if (d < 0.0) d += 1.0;
265 j = d*nazi + 0.5;
266 if (j >= nazi) j = 0;
267 return(i*nazi + j);
268 }
269
270
271 int
272 my_default( /* default illum action */
273 OBJREC *ob,
274 struct illum_args *il,
275 char *nm
276 )
277 {
278 sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
279 nm, ofun[ob->otype].funame, ob->oname);
280 error(WARNING, errmsg);
281 printobj(il->altmat, ob);
282 return(1);
283 }
284
285
286 int
287 my_face( /* make an illum face */
288 OBJREC *ob,
289 struct illum_args *il,
290 char *nm
291 )
292 {
293 int dim[2];
294 int n, nalt, nazi, alti;
295 double sp[2], r1, r2;
296 int h;
297 FVECT dn, org, dir;
298 FVECT u, v;
299 double ur[2], vr[2];
300 MAT4 xfm;
301 int nallow;
302 FACE *fa;
303 int i, j;
304 /* get/check arguments */
305 fa = getface(ob);
306 if (fa->area == 0.0) {
307 freeface(ob);
308 return(my_default(ob, il, nm));
309 }
310 /* set up sampling */
311 if (il->sd != NULL) {
312 if (!getBSDF_xfm(xfm, fa->norm, il->udir)) {
313 objerror(ob, WARNING, "illegal up direction");
314 freeface(ob);
315 return(my_default(ob, il, nm));
316 }
317 n = il->sd->ninc;
318 } else {
319 if (il->sampdens <= 0) {
320 nalt = nazi = 1; /* diffuse assumption */
321 } else {
322 n = PI * il->sampdens;
323 nalt = sqrt(n/PI) + .5;
324 nazi = PI*nalt + .5;
325 }
326 n = nazi*nalt;
327 }
328 newdist(n);
329 /* take first edge >= sqrt(area) */
330 for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
331 u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
332 u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
333 u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
334 if ((r1 = DOT(u,u)) >= fa->area-FTINY)
335 break;
336 }
337 if (i < fa->nv) { /* got one! -- let's align our axes */
338 r2 = 1.0/sqrt(r1);
339 u[0] *= r2; u[1] *= r2; u[2] *= r2;
340 fcross(v, fa->norm, u);
341 } else /* oh well, we'll just have to wing it */
342 mkaxes(u, v, fa->norm);
343 /* now, find limits in (u,v) coordinates */
344 ur[0] = vr[0] = FHUGE;
345 ur[1] = vr[1] = -FHUGE;
346 for (i = 0; i < fa->nv; i++) {
347 r1 = DOT(VERTEX(fa,i),u);
348 if (r1 < ur[0]) ur[0] = r1;
349 if (r1 > ur[1]) ur[1] = r1;
350 r2 = DOT(VERTEX(fa,i),v);
351 if (r2 < vr[0]) vr[0] = r2;
352 if (r2 > vr[1]) vr[1] = r2;
353 }
354 dim[0] = random();
355 /* sample polygon */
356 nallow = 5*n*il->nsamps;
357 for (dim[1] = 0; dim[1] < n; dim[1]++)
358 for (i = 0; i < il->nsamps; i++) {
359 /* randomize direction */
360 h = ilhash(dim, 2) + i;
361 if (il->sd != NULL) {
362 r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
363 } else {
364 multisamp(sp, 2, urand(h));
365 alti = dim[1]/nazi;
366 r1 = (alti + sp[0])/nalt;
367 r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
368 flatdir(dn, r1, r2);
369 for (j = 0; j < 3; j++)
370 dir[j] = -dn[0]*u[j] - dn[1]*v[j] -
371 dn[2]*fa->norm[j];
372 }
373 /* randomize location */
374 do {
375 multisamp(sp, 2, urand(h+4862+nallow));
376 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
377 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
378 for (j = 0; j < 3; j++)
379 org[j] = r1*u[j] + r2*v[j]
380 + fa->offset*fa->norm[j];
381 } while (!inface(org, fa) && nallow-- > 0);
382 if (nallow < 0) {
383 objerror(ob, WARNING, "bad aspect");
384 rayclean();
385 freeface(ob);
386 return(my_default(ob, il, nm));
387 }
388 if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY)
389 r1 = -1.0001*il->thick - 5.*FTINY;
390 else
391 r1 = 5.*FTINY;
392 for (j = 0; j < 3; j++)
393 org[j] += r1*fa->norm[j];
394 /* send sample */
395 raysamp(dim[1], org, dir);
396 }
397 /* add in direct component? */
398 if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
399 MAT4 ixfm;
400 if (il->sd == NULL) {
401 for (i = 3; i--; ) {
402 ixfm[i][0] = u[i];
403 ixfm[i][1] = v[i];
404 ixfm[i][2] = fa->norm[i];
405 ixfm[i][3] = 0.;
406 }
407 ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
408 ixfm[3][3] = 1.;
409 } else {
410 if (!invmat4(ixfm, xfm))
411 objerror(ob, INTERNAL,
412 "cannot invert BSDF transform");
413 if (!(il->flags & IL_LIGHT))
414 new_discount();
415 }
416 dim[0] = random();
417 nallow = 10*il->nsamps;
418 for (i = 0; i < il->nsamps; i++) {
419 /* randomize location */
420 h = dim[0] + samplendx++;
421 do {
422 multisamp(sp, 2, urand(h+nallow));
423 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
424 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
425 for (j = 0; j < 3; j++)
426 org[j] = r1*u[j] + r2*v[j]
427 + fa->offset*fa->norm[j];
428 } while (!inface(org, fa) && nallow-- > 0);
429 if (nallow < 0) {
430 objerror(ob, WARNING, "bad aspect");
431 rayclean();
432 freeface(ob);
433 return(my_default(ob, il, nm));
434 }
435 /* sample source rays */
436 srcsamps(il, org, fa->norm, ixfm);
437 }
438 }
439 /* wait for all rays to finish */
440 rayclean();
441 if (il->sd != NULL) { /* run distribution through BSDF */
442 nalt = sqrt(il->sd->nout/PI) + .5;
443 nazi = PI*nalt + .5;
444 redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm);
445 done_discount();
446 if (!il->sampdens)
447 il->sampdens = nalt*nazi/PI + .999;
448 }
449 /* write out the face and its distribution */
450 if (average(il, distarr, n)) {
451 if (il->sampdens > 0)
452 flatout(il, distarr, nalt, nazi, u, v, fa->norm);
453 illumout(il, ob);
454 } else
455 printobj(il->altmat, ob);
456 /* clean up */
457 freeface(ob);
458 return(0);
459 }
460
461
462 int
463 my_sphere( /* make an illum sphere */
464 register OBJREC *ob,
465 struct illum_args *il,
466 char *nm
467 )
468 {
469 int dim[3];
470 int n, nalt, nazi;
471 double sp[4], r1, r2, r3;
472 FVECT org, dir;
473 FVECT u, v;
474 register int i, j;
475 /* check arguments */
476 if (ob->oargs.nfargs != 4)
477 objerror(ob, USER, "bad # of arguments");
478 /* set up sampling */
479 if (il->sampdens <= 0)
480 nalt = nazi = 1;
481 else {
482 n = 4.*PI * il->sampdens;
483 nalt = sqrt(2./PI*n) + .5;
484 nazi = PI/2.*nalt + .5;
485 }
486 if (il->sd != NULL)
487 objerror(ob, WARNING, "BSDF ignored");
488 n = nalt*nazi;
489 newdist(n);
490 dim[0] = random();
491 /* sample sphere */
492 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
493 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
494 for (i = 0; i < il->nsamps; i++) {
495 /* next sample point */
496 multisamp(sp, 4, urand(ilhash(dim,3)+i));
497 /* random direction */
498 r1 = (dim[1] + sp[0])/nalt;
499 r2 = (dim[2] + sp[1] - .5)/nazi;
500 rounddir(dir, r1, r2);
501 /* random location */
502 mkaxes(u, v, dir); /* yuck! */
503 r3 = sqrt(sp[2]);
504 r2 = 2.*PI*sp[3];
505 r1 = r3*ob->oargs.farg[3]*cos(r2);
506 r2 = r3*ob->oargs.farg[3]*sin(r2);
507 r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
508 for (j = 0; j < 3; j++) {
509 org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
510 r3*dir[j];
511 dir[j] = -dir[j];
512 }
513 /* send sample */
514 raysamp(dim[1]*nazi+dim[2], org, dir);
515 }
516 /* wait for all rays to finish */
517 rayclean();
518 /* write out the sphere and its distribution */
519 if (average(il, distarr, n)) {
520 if (il->sampdens > 0)
521 roundout(il, distarr, nalt, nazi);
522 else
523 objerror(ob, WARNING, "diffuse distribution");
524 illumout(il, ob);
525 } else
526 printobj(il->altmat, ob);
527 /* clean up */
528 return(1);
529 }
530
531
532 int
533 my_ring( /* make an illum ring */
534 OBJREC *ob,
535 struct illum_args *il,
536 char *nm
537 )
538 {
539 int dim[2];
540 int n, nalt, nazi, alti;
541 double sp[2], r1, r2, r3;
542 int h;
543 FVECT dn, org, dir;
544 FVECT u, v;
545 MAT4 xfm;
546 CONE *co;
547 int i, j;
548 /* get/check arguments */
549 co = getcone(ob, 0);
550 /* set up sampling */
551 if (il->sd != NULL) {
552 if (!getBSDF_xfm(xfm, co->ad, il->udir)) {
553 objerror(ob, WARNING, "illegal up direction");
554 freecone(ob);
555 return(my_default(ob, il, nm));
556 }
557 n = il->sd->ninc;
558 } else {
559 if (il->sampdens <= 0) {
560 nalt = nazi = 1; /* diffuse assumption */
561 } else {
562 n = PI * il->sampdens;
563 nalt = sqrt(n/PI) + .5;
564 nazi = PI*nalt + .5;
565 }
566 n = nazi*nalt;
567 }
568 newdist(n);
569 mkaxes(u, v, co->ad);
570 dim[0] = random();
571 /* sample disk */
572 for (dim[1] = 0; dim[1] < n; dim[1]++)
573 for (i = 0; i < il->nsamps; i++) {
574 /* next sample point */
575 h = ilhash(dim,2) + i;
576 /* randomize direction */
577 if (il->sd != NULL) {
578 r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
579 } else {
580 multisamp(sp, 2, urand(h));
581 alti = dim[1]/nazi;
582 r1 = (alti + sp[0])/nalt;
583 r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
584 flatdir(dn, r1, r2);
585 for (j = 0; j < 3; j++)
586 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
587 }
588 /* randomize location */
589 multisamp(sp, 2, urand(h+8371));
590 r3 = sqrt(CO_R0(co)*CO_R0(co) +
591 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
592 r2 = 2.*PI*sp[1];
593 r1 = r3*cos(r2);
594 r2 = r3*sin(r2);
595 if (il->sd != NULL && DOT(dir, co->ad) < -FTINY)
596 r3 = -1.0001*il->thick - 5.*FTINY;
597 else
598 r3 = 5.*FTINY;
599 for (j = 0; j < 3; j++)
600 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
601 r3*co->ad[j];
602 /* send sample */
603 raysamp(dim[1], org, dir);
604 }
605 /* add in direct component? */
606 if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
607 MAT4 ixfm;
608 if (il->sd == NULL) {
609 for (i = 3; i--; ) {
610 ixfm[i][0] = u[i];
611 ixfm[i][1] = v[i];
612 ixfm[i][2] = co->ad[i];
613 ixfm[i][3] = 0.;
614 }
615 ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
616 ixfm[3][3] = 1.;
617 } else {
618 if (!invmat4(ixfm, xfm))
619 objerror(ob, INTERNAL,
620 "cannot invert BSDF transform");
621 if (!(il->flags & IL_LIGHT))
622 new_discount();
623 }
624 dim[0] = random();
625 for (i = 0; i < il->nsamps; i++) {
626 /* randomize location */
627 h = dim[0] + samplendx++;
628 multisamp(sp, 2, urand(h));
629 r3 = sqrt(CO_R0(co)*CO_R0(co) +
630 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
631 r2 = 2.*PI*sp[1];
632 r1 = r3*cos(r2);
633 r2 = r3*sin(r2);
634 for (j = 0; j < 3; j++)
635 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j];
636 /* sample source rays */
637 srcsamps(il, org, co->ad, ixfm);
638 }
639 }
640 /* wait for all rays to finish */
641 rayclean();
642 if (il->sd != NULL) { /* run distribution through BSDF */
643 nalt = sqrt(il->sd->nout/PI) + .5;
644 nazi = PI*nalt + .5;
645 redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
646 done_discount();
647 if (!il->sampdens)
648 il->sampdens = nalt*nazi/PI + .999;
649 }
650 /* write out the ring and its distribution */
651 if (average(il, distarr, n)) {
652 if (il->sampdens > 0)
653 flatout(il, distarr, nalt, nazi, u, v, co->ad);
654 illumout(il, ob);
655 } else
656 printobj(il->altmat, ob);
657 /* clean up */
658 freecone(ob);
659 return(1);
660 }
661
662
663 void
664 redistribute( /* pass distarr ray sums through BSDF */
665 struct BSDF_data *b,
666 int nalt,
667 int nazi,
668 FVECT u,
669 FVECT v,
670 FVECT w,
671 MAT4 xm
672 )
673 {
674 int nout = 0;
675 MAT4 mymat, inmat;
676 COLORV *idist;
677 COLORV *cp;
678 FVECT dv;
679 double wt;
680 int i, j, k, c, o;
681 COLOR col, cinc;
682 /* copy incoming distribution */
683 if (b->ninc > distsiz)
684 error(INTERNAL, "error 1 in redistribute");
685 idist = (COLORV *)malloc(sizeof(COLOR)*b->ninc);
686 if (idist == NULL)
687 error(SYSTEM, "out of memory in redistribute");
688 memcpy(idist, distarr, sizeof(COLOR)*b->ninc);
689 /* compose direction transform */
690 for (i = 3; i--; ) {
691 mymat[i][0] = u[i];
692 mymat[i][1] = v[i];
693 mymat[i][2] = w[i];
694 mymat[i][3] = 0.;
695 }
696 mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
697 mymat[3][3] = 1.;
698 if (xm != NULL)
699 multmat4(mymat, xm, mymat);
700 for (i = 3; i--; ) { /* make sure it's normalized */
701 wt = 1./sqrt( mymat[0][i]*mymat[0][i] +
702 mymat[1][i]*mymat[1][i] +
703 mymat[2][i]*mymat[2][i] );
704 for (j = 3; j--; )
705 mymat[j][i] *= wt;
706 }
707 if (!invmat4(inmat, mymat)) /* need inverse as well */
708 error(INTERNAL, "cannot invert BSDF transform");
709 newdist(nalt*nazi); /* resample distribution */
710 for (i = b->ninc; i--; ) {
711 int direct_out = -1;
712 COLOR cdir;
713 getBSDF_incvec(dv, b, i); /* compute incident irrad. */
714 multv3(dv, dv, mymat);
715 if (dv[2] < 0.0) {
716 dv[0] = -dv[0]; dv[1] = -dv[1]; dv[2] = -dv[2];
717 direct_out += (direct_discount != NULL);
718 }
719 wt = getBSDF_incohm(b, i);
720 wt *= dv[2]; /* solid_angle*cosine(theta) */
721 cp = &idist[3*i];
722 copycolor(cinc, cp);
723 scalecolor(cinc, wt);
724 if (!direct_out) { /* discount direct contr. */
725 cp = &direct_discount[3*i];
726 copycolor(cdir, cp);
727 scalecolor(cdir, -wt);
728 if (b->nout != b->ninc)
729 direct_out = flatindex(dv, nalt, nazi);
730 else
731 direct_out = i; /* assumes dist. mirroring */
732 }
733 for (k = nalt; k--; ) /* loop over distribution */
734 for (j = nazi; j--; ) {
735 int rstart = random();
736 for (c = NBSDFSAMPS; c--; ) {
737 double sp[2];
738 multisamp(sp, 2, urand(rstart+c));
739 flatdir(dv, (k + sp[0])/nalt,
740 (j + .5 - sp[1])/nazi);
741 multv3(dv, dv, inmat);
742 /* evaluate BSDF @ outgoing */
743 o = getBSDF_outndx(b, dv);
744 if (o < 0) {
745 nout++;
746 continue;
747 }
748 wt = BSDF_value(b, i, o) * (1./NBSDFSAMPS);
749 copycolor(col, cinc);
750 if (b->nout != b->ninc)
751 o = k*nazi + j;
752 if (o == direct_out)
753 addcolor(col, cdir); /* minus direct */
754 scalecolor(col, wt);
755 cp = &distarr[3*(k*nazi + j)];
756 addcolor(cp, col); /* sum into distribution */
757 }
758 }
759 }
760 free(idist); /* free temp space */
761 if (nout) {
762 sprintf(errmsg, "missing %.1f%% of BSDF directions",
763 100.*nout/(b->ninc*nalt*nazi*NBSDFSAMPS));
764 error(WARNING, errmsg);
765 }
766 }