ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapmat.c
Revision: 2.2
Committed: Wed Apr 22 15:50:44 2015 UTC (9 years, 1 month ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.1: +375 -138 lines
Log Message:
bsdfPhotonScatter() now uses SDsampBSDF() from BSDF library rather than
homebrew sampling code...

File Contents

# Content
1 /*
2 ==================================================================
3 Photon map support routines for scattering by materials.
4
5 Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
6 (c) Fraunhofer Institute for Solar Energy Systems,
7 Lucerne University of Applied Sciences & Arts
8 ==================================================================
9
10 $Id: pmapmat.c,v 4.14 2015/04/22 15:48:42 taschreg Exp taschreg $
11 */
12
13
14
15 #include "pmapmat.h"
16 #include "pmapdata.h"
17 #include "pmaprand.h"
18 #include "otypes.h"
19 #include "data.h"
20 #include "func.h"
21 #include "bsdf.h"
22 #include <math.h>
23
24
25
26 /* Stuff ripped off from material modules */
27 #define MAXITER 10
28 #define SP_REFL 01
29 #define SP_TRAN 02
30 #define SP_PURE 04
31 #define SP_FLAT 010
32 #define SP_BADU 040
33 #define MLAMBDA 500
34 #define RINDEX 1.52
35 #define FRESNE(ci) (exp(-5.85*(ci)) - 0.00287989916)
36
37
38
39 typedef struct {
40 OBJREC *mp;
41 RAY *rp;
42 short specfl;
43 COLOR mcolor, scolor;
44 FVECT vrefl, prdir, pnorm;
45 double alpha2, rdiff, rspec, trans, tdiff, tspec, pdot;
46 } NORMDAT;
47
48 typedef struct {
49 OBJREC *mp;
50 RAY *rp;
51 short specfl;
52 COLOR mcolor, scolor;
53 FVECT vrefl, prdir, u, v, pnorm;
54 double u_alpha, v_alpha, rdiff, rspec, trans, tdiff, tspec, pdot;
55 } ANISODAT;
56
57 typedef struct {
58 OBJREC *mp;
59 RAY *pr;
60 FVECT pnorm;
61 FVECT vray;
62 double sr_vpsa [2];
63 RREAL toloc [3][3];
64 RREAL fromloc [3][3];
65 double thick;
66 SDData *sd;
67 COLOR runsamp;
68 COLOR rdiff;
69 COLOR tunsamp;
70 COLOR tdiff;
71 } BSDFDAT;
72
73
74
75 extern const SDCDst SDemptyCD;
76
77 /* Per-material scattering function dispatch table; return value is usually
78 * zero, indicating photon termination */
79 int (*photonScatter [NUMOTYPE]) (OBJREC*, RAY*);
80
81 /* List of antimatter sensor modifier names and associated object set */
82 char *photonSensorList [MAXSET + 1] = {NULL};
83 static OBJECT photonSensorSet [MAXSET + 1] = {0};
84
85
86
87 /* ================ General support routines ================ */
88
89
90 void photonRay (const RAY *rayIn, RAY *rayOut,
91 int rayOutType, COLOR fluxAtten)
92 /* Spawn a new photon ray from a previous one; this is effectively a
93 * customised rayorigin().
94 * A SPECULAR rayOutType flags this photon as _caustic_ for subsequent hits.
95 * It is preserved for transferred rays (of type PMAP_XFER).
96 * fluxAtten specifies the RGB attenuation of the photon flux effected by
97 * the scattering material. The outgoing flux is then normalised to maintain
98 * a uniform average of 1 over RGB. If fluxAtten == NULL, the flux remains
99 * unchanged for the outgoing photon. fluxAtten is ignored for transferred
100 * rays.
101 * The ray direction is preserved for transferred rays, and undefined for
102 * scattered rays and must be subsequently set by the caller. */
103 {
104 rayorigin(rayOut, rayOutType, rayIn, NULL);
105
106 /* Transfer flux */
107 copycolor(rayOut -> rcol, rayIn -> rcol);
108
109 /* Copy caustic flag & direction for transferred rays */
110 if (rayOutType == PMAP_XFER) {
111 /* rayOut -> rtype |= rayIn -> rtype & SPECULAR; */
112 rayOut -> rtype |= rayIn -> rtype;
113 VCOPY(rayOut -> rdir, rayIn -> rdir);
114 }
115 else if (fluxAtten) {
116 /* Attenuate and normalise flux for scattered rays */
117 multcolor(rayOut -> rcol, fluxAtten);
118 colorNorm(rayOut -> rcol);
119 }
120
121 /* Propagate index of emitting light source */
122 rayOut -> rsrc = rayIn -> rsrc;
123 }
124
125
126
127 static void addPhotons (const RAY *r)
128 /* Insert photon hits, where applicable */
129 {
130 if (!r -> rlvl)
131 /* Add direct photon map at primary hitpoint */
132 addPhoton(directPmap, r);
133 else {
134 /* Add global or precomputed photon map at indirect hitpoint */
135 addPhoton(preCompPmap ? preCompPmap : globalPmap, r);
136
137 /* Store caustic photon if specular flag set */
138 if (PMAP_CAUSTICRAY(r))
139 addPhoton(causticPmap, r);
140
141 /* Store in contribution photon map */
142 addPhoton(contribPmap, r);
143 }
144 }
145
146
147
148 void getPhotonSensors (char **sensorList)
149 /* Find antimatter geometry declared as photon sensors */
150 {
151 OBJECT i;
152 OBJREC *obj;
153 char **lp;
154
155 /* Init sensor set */
156 photonSensorSet [0] = 0;
157
158 if (!sensorList [0])
159 return;
160
161 for (i = 0; i < nobjects; i++) {
162 obj = objptr(i);
163
164 /* Insert object in sensor set if it's in the specified sensor list
165 * and of type antimatter */
166 for (lp = sensorList; *lp; lp++) {
167 if (!strcmp(obj -> oname, *lp)) {
168 if (obj -> otype != MAT_CLIP) {
169 sprintf(errmsg, "photon sensor modifier %s is not antimatter",
170 obj -> oname);
171 error(USER, errmsg);
172 }
173
174 if (photonSensorSet [0] >= AMBLLEN)
175 error(USER, "too many photon sensor modifiers");
176
177 insertelem(photonSensorSet, i);
178 }
179 }
180 }
181
182 if (!photonSensorSet [0])
183 error(USER, "no photon sensors found");
184 }
185
186
187
188 /* ================ Material specific scattering routines ================ */
189
190
191 static int isoSpecPhotonScatter (NORMDAT *nd, RAY *rayOut)
192 /* Generate direction for isotropically specularly reflected
193 or transmitted ray. Returns 1 if successful. */
194 {
195 FVECT u, v, h;
196 RAY *rayIn = nd -> rp;
197 double d, d2, sinp, cosp;
198 int niter, i = 0;
199
200 /* Set up sample coordinates */
201 do {
202 v [0] = v [1] = v [2] = 0;
203 v [i++] = 1;
204 fcross(u, v, nd -> pnorm);
205 } while (normalize(u) < FTINY);
206
207 fcross(v, nd -> pnorm, u);
208
209 if (nd -> specfl & SP_REFL) {
210 /* Specular reflection; make MAXITER attempts at getting a ray */
211
212 for (niter = 0; niter < MAXITER; niter++) {
213 d = 2 * PI * pmapRandom(scatterState);
214 cosp = cos(d);
215 sinp = sin(d);
216 d2 = pmapRandom(scatterState);
217 d = d2 <= FTINY ? 1 : sqrt(nd -> alpha2 * -log(d2));
218
219 for (i = 0; i < 3; i++)
220 h [i] = nd -> pnorm [i] + d * (cosp * u [i] + sinp * v [i]);
221
222 d = -2 * DOT(h, rayIn -> rdir) / (1 + d * d);
223 VSUM(rayOut -> rdir, rayIn -> rdir, h, d);
224
225 if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY)
226 return 1;
227 }
228
229 return 0;
230 }
231
232 else {
233 /* Specular transmission; make MAXITER attempts at getting a ray */
234
235 for (niter = 0; niter < MAXITER; niter++) {
236 d = 2 * PI * pmapRandom(scatterState);
237 cosp = cos(d);
238 sinp = sin(d);
239 d2 = pmapRandom(scatterState);
240 d = d2 <= FTINY ? 1 : sqrt(-log(d2) * nd -> alpha2);
241
242 for (i = 0; i < 3; i++)
243 rayOut -> rdir [i] = nd -> prdir [i] +
244 d * (cosp * u [i] + sinp * v [i]);
245
246 if (DOT(rayOut -> rdir, rayIn -> ron) < -FTINY) {
247 normalize(rayOut -> rdir);
248 return 1;
249 }
250 }
251
252 return 0;
253 }
254 }
255
256
257
258 static void diffPhotonScatter (FVECT normal, RAY* rayOut)
259 /* Generate cosine-weighted direction for diffuse ray */
260 {
261 const RREAL cosThetaSqr = pmapRandom(scatterState),
262 cosTheta = sqrt(cosThetaSqr),
263 sinTheta = sqrt(1 - cosThetaSqr),
264 phi = 2 * PI * pmapRandom(scatterState),
265 du = cos(phi) * sinTheta, dv = sin(phi) * sinTheta;
266 FVECT u, v;
267 int i = 0;
268
269 /* Set up sample coordinates */
270 do {
271 v [0] = v [1] = v [2] = 0;
272 v [i++] = 1;
273 fcross(u, v, normal);
274 } while (normalize(u) < FTINY);
275
276 fcross(v, normal, u);
277
278 /* Convert theta & phi to cartesian */
279 for (i = 0; i < 3; i++)
280 rayOut -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * normal [i];
281
282 normalize(rayOut -> rdir);
283 }
284
285
286
287 static int normalPhotonScatter (OBJREC *mat, RAY *rayIn)
288 /* Generate new photon ray for isotropic material and recurse */
289 {
290 NORMDAT nd;
291 int i, hastexture;
292 float xi, albedo, prdiff, ptdiff, prspec, ptspec;
293 double d, fresnel;
294 RAY rayOut;
295
296 if (mat -> oargs.nfargs != (mat -> otype == MAT_TRANS ? 7 : 5))
297 objerror(mat, USER, "bad number of arguments");
298
299 /* Check for back side; reorient if back is visible */
300 if (rayIn -> rod < 0)
301 if (!backvis && mat -> otype != MAT_TRANS)
302 return 0;
303 else {
304 /* Get modifiers */
305 raytexture(rayIn, mat -> omod);
306 flipsurface(rayIn);
307 }
308 else raytexture(rayIn, mat -> omod);
309
310 nd.rp = rayIn;
311
312 /* Get material color */
313 copycolor(nd.mcolor, mat -> oargs.farg);
314
315 /* Get roughness */
316 nd.specfl = 0;
317 nd.alpha2 = mat -> oargs.farg [4];
318
319 if ((nd.alpha2 *= nd.alpha2) <= FTINY)
320 nd.specfl |= SP_PURE;
321
322 if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype))
323 nd.specfl |= SP_FLAT;
324
325 /* Perturb normal */
326 if ((hastexture = DOT(rayIn -> pert, rayIn -> pert)) > sqr(FTINY))
327 nd.pdot = raynormal(nd.pnorm, rayIn);
328 else {
329 VCOPY(nd.pnorm, rayIn -> ron);
330 nd.pdot = rayIn -> rod;
331 }
332
333 nd.pdot = max(nd.pdot, .001);
334
335 /* Modify material color */
336 multcolor(nd.mcolor, rayIn -> pcol);
337 nd.rspec = mat -> oargs.farg [3];
338
339 /* Approximate Fresnel term */
340 if (nd.specfl & SP_PURE && nd.rspec > FTINY) {
341 fresnel = FRESNE(rayIn -> rod);
342 nd.rspec += fresnel * (1 - nd.rspec);
343 }
344 else fresnel = 0;
345
346 /* Transmission params */
347 if (mat -> otype == MAT_TRANS) {
348 nd.trans = mat -> oargs.farg [5] * (1 - nd.rspec);
349 nd.tspec = nd.trans * mat -> oargs.farg [6];
350 nd.tdiff = nd.trans - nd.tspec;
351 }
352 else nd.tdiff = nd.tspec = nd.trans = 0;
353
354 /* Specular reflection params */
355 if (nd.rspec > FTINY) {
356 /* Specular color */
357 if (mat -> otype != MAT_METAL)
358 setcolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec);
359 else if (fresnel > FTINY) {
360 d = nd.rspec * (1 - fresnel);
361 for (i = 0; i < 3; i++)
362 nd.scolor [i] = fresnel + nd.mcolor [i] * d;
363 }
364 else {
365 copycolor(nd.scolor, nd.mcolor);
366 scalecolor(nd.scolor, nd.rspec);
367 }
368 }
369 else setcolor(nd.scolor, 0, 0, 0);
370
371 /* Diffuse reflection params */
372 nd.rdiff = 1 - nd.trans - nd.rspec;
373
374 /* Set up probabilities */
375 prdiff = ptdiff = ptspec = colorAvg(nd.mcolor);
376 prdiff *= nd.rdiff;
377 ptdiff *= nd.tdiff;
378 prspec = colorAvg(nd.scolor);
379 ptspec *= nd.tspec;
380 albedo = prdiff + ptdiff + prspec + ptspec;
381
382 /* Insert direct and indirect photon hits if diffuse component */
383 if (prdiff > FTINY || ptdiff > FTINY)
384 addPhotons(rayIn);
385
386 xi = pmapRandom(rouletteState);
387
388 if (xi > albedo)
389 /* Absorbed */
390 return 0;
391
392 if (xi > (albedo -= prspec)) {
393 /* Specular reflection */
394 nd.specfl |= SP_REFL;
395
396 if (nd.specfl & SP_PURE) {
397 /* Perfect specular reflection */
398 for (i = 0; i < 3; i++) {
399 /* Reflected ray */
400 nd.vrefl [i] = rayIn -> rdir [i] + 2 * nd.pdot * nd.pnorm [i];
401 }
402
403 /* Penetration? */
404 if (hastexture && DOT(nd.vrefl, rayIn -> ron) <= FTINY)
405 for (i = 0; i < 3; i++) {
406 /* Safety measure */
407 nd.vrefl [i] = rayIn -> rdir [i] +
408 2 * rayIn -> rod * rayIn -> ron [i];
409 }
410
411 VCOPY(rayOut.rdir, nd.vrefl);
412 }
413
414 else if (!isoSpecPhotonScatter(&nd, &rayOut))
415 return 0;
416
417 photonRay(rayIn, &rayOut, PMAP_SPECREFL, nd.scolor);
418 }
419
420 else if (xi > (albedo -= ptspec)) {
421 /* Specular transmission */
422 nd.specfl |= SP_TRAN;
423
424 if (hastexture) {
425 /* Perturb */
426 for (i = 0; i < 3; i++)
427 nd.prdir [i] = rayIn -> rdir [i] - rayIn -> pert [i];
428
429 if (DOT(nd.prdir, rayIn -> ron) < -FTINY)
430 normalize(nd.prdir);
431 else VCOPY(nd.prdir, rayIn -> rdir);
432 }
433 else VCOPY(nd.prdir, rayIn -> rdir);
434
435 if ((nd.specfl & (SP_TRAN | SP_PURE)) == (SP_TRAN | SP_PURE))
436 /* Perfect specular transmission */
437 VCOPY(rayOut.rdir, nd.prdir);
438 else if (!isoSpecPhotonScatter(&nd, &rayOut))
439 return 0;
440
441 photonRay(rayIn, &rayOut, PMAP_SPECTRANS, nd.mcolor);
442 }
443
444 else if (xi > (albedo -= prdiff)) {
445 /* Diffuse reflection */
446 photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor);
447 diffPhotonScatter(hastexture ? nd.pnorm : rayIn -> ron, &rayOut);
448 }
449
450 else {
451 /* Diffuse transmission */
452 flipsurface(rayIn);
453 photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor);
454
455 if (hastexture) {
456 FVECT bnorm;
457 bnorm [0] = -nd.pnorm [0];
458 bnorm [1] = -nd.pnorm [1];
459 bnorm [2] = -nd.pnorm [2];
460 diffPhotonScatter(bnorm, &rayOut);
461 }
462 else diffPhotonScatter(rayIn -> ron, &rayOut);
463 }
464
465 tracePhoton(&rayOut);
466 return 0;
467 }
468
469
470
471 static void getacoords (ANISODAT *np)
472 /* Set up coordinate system for anisotropic sampling; cloned from aniso.c */
473 {
474 MFUNC *mf;
475 int i;
476
477 mf = getfunc(np->mp, 3, 0x7, 1);
478 setfunc(np->mp, np->rp);
479 errno = 0;
480
481 for (i = 0; i < 3; i++)
482 np->u[i] = evalue(mf->ep[i]);
483
484 if ((errno == EDOM) | (errno == ERANGE)) {
485 objerror(np->mp, WARNING, "compute error");
486 np->specfl |= SP_BADU;
487 return;
488 }
489
490 if (mf->fxp != &unitxf)
491 multv3(np->u, np->u, mf->fxp->xfm);
492
493 fcross(np->v, np->pnorm, np->u);
494
495 if (normalize(np->v) == 0.0) {
496 objerror(np->mp, WARNING, "illegal orientation vector");
497 np->specfl |= SP_BADU;
498 return;
499 }
500
501 fcross(np->u, np->v, np->pnorm);
502 }
503
504
505
506 static int anisoSpecPhotonScatter (ANISODAT *nd, RAY *rayOut)
507 /* Generate direction for anisotropically specularly reflected
508 or transmitted ray. Returns 1 if successful. */
509 {
510 FVECT h;
511 double d, d2, sinp, cosp;
512 int niter, i;
513 RAY *rayIn = nd -> rp;
514
515 if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype))
516 nd -> specfl |= SP_FLAT;
517
518 /* set up coordinates */
519 getacoords(nd);
520
521 if (rayOut -> rtype & TRANS) {
522 /* Specular transmission */
523
524 if (DOT(rayIn -> pert, rayIn -> pert) <= FTINY * FTINY)
525 VCOPY(nd -> prdir, rayIn -> rdir);
526 else {
527 /* perturb */
528 for (i = 0; i < 3; i++)
529 nd -> prdir [i] = rayIn -> rdir [i] - rayIn -> pert [i];
530
531 if (DOT(nd -> prdir, rayIn -> ron) < -FTINY)
532 normalize(nd -> prdir);
533 else VCOPY(nd -> prdir, rayIn -> rdir);
534 }
535
536 /* Make MAXITER attempts at getting a ray */
537 for (niter = 0; niter < MAXITER; niter++) {
538 d = 2 * PI * pmapRandom(scatterState);
539 cosp = cos(d) * nd -> u_alpha;
540 sinp = sin(d) * nd -> v_alpha;
541 d = sqrt(sqr(cosp) + sqr(sinp));
542 cosp /= d;
543 sinp /= d;
544 d2 = pmapRandom(scatterState);
545 d = d2 <= FTINY ? 1
546 : sqrt(-log(d2) /
547 (sqr(cosp) / sqr(nd -> u_alpha) +
548 sqr(sinp) / (nd -> v_alpha * nd -> u_alpha)));
549
550 for (i = 0; i < 3; i++)
551 rayOut -> rdir [i] = nd -> prdir [i] + d *
552 (cosp * nd -> u [i] + sinp * nd -> v [i]);
553
554 if (DOT(rayOut -> rdir, rayIn -> ron) < -FTINY) {
555 normalize(rayOut -> rdir);
556 return 1;
557 }
558 }
559
560 return 0;
561 }
562
563 else {
564 /* Specular reflection */
565
566 /* Make MAXITER attempts at getting a ray */
567 for (niter = 0; niter < MAXITER; niter++) {
568 d = 2 * PI * pmapRandom(scatterState);
569 cosp = cos(d) * nd -> u_alpha;
570 sinp = sin(d) * nd -> v_alpha;
571 d = sqrt(sqr(cosp) + sqr(sinp));
572 cosp /= d;
573 sinp /= d;
574 d2 = pmapRandom(scatterState);
575 d = d2 <= FTINY ? 1
576 : sqrt(-log(d2) /
577 (sqr(cosp) / sqr(nd -> u_alpha) +
578 sqr(sinp) / (nd -> v_alpha * nd -> v_alpha)));
579
580 for (i = 0; i < 3; i++)
581 h [i] = nd -> pnorm [i] +
582 d * (cosp * nd -> u [i] + sinp * nd -> v [i]);
583
584 d = -2 * DOT(h, rayIn -> rdir) / (1 + d * d);
585 VSUM(rayOut -> rdir, rayIn -> rdir, h, d);
586
587 if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY)
588 return 1;
589 }
590
591 return 0;
592 }
593 }
594
595
596
597 static int anisoPhotonScatter (OBJREC *mat, RAY *rayIn)
598 /* Generate new photon ray for anisotropic material and recurse */
599 {
600 ANISODAT nd;
601 float xi, albedo, prdiff, ptdiff, prspec, ptspec;
602 RAY rayOut;
603
604 if (mat -> oargs.nfargs != (mat -> otype == MAT_TRANS2 ? 8 : 6))
605 objerror(mat, USER, "bad number of real arguments");
606
607 nd.rp = rayIn;
608 nd.mp = objptr(rayIn -> ro -> omod);
609
610 /* get material color */
611 copycolor(nd.mcolor, mat -> oargs.farg);
612
613 /* get roughness */
614 nd.specfl = 0;
615 nd.u_alpha = mat -> oargs.farg [4];
616 nd.v_alpha = mat -> oargs.farg [5];
617 if (nd.u_alpha < FTINY || nd.v_alpha <= FTINY)
618 objerror(mat, USER, "roughness too small");
619
620 /* check for back side; reorient if back is visible */
621 if (rayIn -> rod < 0)
622 if (!backvis && mat -> otype != MAT_TRANS2)
623 return 0;
624 else {
625 /* get modifiers */
626 raytexture(rayIn, mat -> omod);
627 flipsurface(rayIn);
628 }
629 else raytexture(rayIn, mat -> omod);
630
631 /* perturb normal */
632 nd.pdot = max(raynormal(nd.pnorm, rayIn), .001);
633
634 /* modify material color */
635 multcolor(nd.mcolor, rayIn -> pcol);
636 nd.rspec = mat -> oargs.farg [3];
637
638 /* transmission params */
639 if (mat -> otype == MAT_TRANS2) {
640 nd.trans = mat -> oargs.farg [6] * (1 - nd.rspec);
641 nd.tspec = nd.trans * mat -> oargs.farg [7];
642 nd.tdiff = nd.trans - nd.tspec;
643 if (nd.tspec > FTINY)
644 nd.specfl |= SP_TRAN;
645 }
646 else nd.tdiff = nd.tspec = nd.trans = 0;
647
648 /* specular reflection params */
649 if (nd.rspec > FTINY) {
650 nd.specfl |= SP_REFL;
651
652 /* comput e specular color */
653 if (mat -> otype == MAT_METAL2)
654 copycolor(nd.scolor, nd.mcolor);
655 else setcolor(nd.scolor, 1, 1, 1);
656
657 scalecolor(nd.scolor, nd.rspec);
658 }
659 else setcolor(nd.scolor, 0, 0, 0);
660
661 /* diffuse reflection params */
662 nd.rdiff = 1 - nd.trans - nd.rspec;
663
664 /* Set up probabilities */
665 prdiff = ptdiff = ptspec = colorAvg(nd.mcolor);
666 prdiff *= nd.rdiff;
667 ptdiff *= nd.tdiff;
668 prspec = colorAvg(nd.scolor);
669 ptspec *= nd.tspec;
670 albedo = prdiff + ptdiff + prspec + ptspec;
671
672 /* Insert direct and indirect photon hits if diffuse component */
673 if (prdiff > FTINY || ptdiff > FTINY)
674 addPhotons(rayIn);
675
676 xi = pmapRandom(rouletteState);
677
678 if (xi > albedo)
679 /* Absorbed */
680 return 0;
681
682 if (xi > (albedo -= prspec))
683 /* Specular reflection */
684 if (!(nd.specfl & SP_BADU)) {
685 photonRay(rayIn, &rayOut, PMAP_SPECREFL, nd.scolor);
686
687 if (!anisoSpecPhotonScatter(&nd, &rayOut))
688 return 0;
689 }
690 else return 0;
691
692 else if (xi > (albedo -= ptspec))
693 /* Specular transmission */
694
695 if (!(nd.specfl & SP_BADU)) {
696 /* Specular transmission */
697 photonRay(rayIn, &rayOut, PMAP_SPECTRANS, nd.mcolor);
698
699 if (!anisoSpecPhotonScatter(&nd, &rayOut))
700 return 0;
701 }
702 else return 0;
703
704 else if (xi > (albedo -= prdiff)) {
705 /* Diffuse reflection */
706 photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor);
707 diffPhotonScatter(nd.pnorm, &rayOut);
708 }
709
710 else {
711 /* Diffuse transmission */
712 FVECT bnorm;
713 flipsurface(rayIn);
714 bnorm [0] = -nd.pnorm [0];
715 bnorm [1] = -nd.pnorm [1];
716 bnorm [2] = -nd.pnorm [2];
717
718 photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor);
719 diffPhotonScatter(bnorm, &rayOut);
720 }
721
722 tracePhoton(&rayOut);
723 return 0;
724 }
725
726
727 static double mylog (double x)
728 /* special log for extinction coefficients; cloned from dielectric.c */
729 {
730 if (x < 1e-40)
731 return(-100.);
732
733 if (x >= 1.)
734 return(0.);
735
736 return(log(x));
737 }
738
739
740 static int dielectricPhotonScatter (OBJREC *mat, RAY *rayIn)
741 /* Generate new photon ray for dielectric material and recurse */
742 {
743 double cos1, cos2, nratio, d1, d2, refl;
744 COLOR ctrans, talb;
745 FVECT dnorm;
746 int hastexture, i;
747 RAY rayOut;
748
749 if (mat -> oargs.nfargs != (mat -> otype == MAT_DIELECTRIC ? 5 : 8))
750 objerror(mat, USER, "bad arguments");
751
752 /* get modifiers */
753 raytexture(rayIn, mat -> omod);
754
755 if ((hastexture = DOT(rayIn -> pert, rayIn -> pert)) > FTINY * FTINY)
756 /* Perturb normal */
757 cos1 = raynormal(dnorm, rayIn);
758 else {
759 VCOPY(dnorm, rayIn -> ron);
760 cos1 = rayIn -> rod;
761 }
762
763 /* index of refraction */
764 nratio = mat -> otype ==
765 MAT_DIELECTRIC ? mat -> oargs.farg [3] + mat -> oargs.farg [4] / MLAMBDA
766 : mat -> oargs.farg [3] / mat -> oargs.farg [7];
767
768 if (cos1 < 0) {
769 /* inside */
770 hastexture = -hastexture;
771 cos1 = -cos1;
772 dnorm [0] = -dnorm [0];
773 dnorm [1] = -dnorm [1];
774 dnorm [2] = -dnorm [2];
775 setcolor(rayIn -> cext,
776 -mylog(mat -> oargs.farg [0] * rayIn -> pcol [0]),
777 -mylog(mat -> oargs.farg [1] * rayIn -> pcol [1]),
778 -mylog(mat -> oargs.farg [2] * rayIn -> pcol [2]));
779 setcolor(rayIn -> albedo, 0, 0, 0);
780 rayIn -> gecc = 0;
781
782 if (mat -> otype == MAT_INTERFACE) {
783 setcolor(ctrans,
784 -mylog(mat -> oargs.farg [4] * rayIn -> pcol [0]),
785 -mylog(mat -> oargs.farg [5] * rayIn -> pcol [1]),
786 -mylog(mat -> oargs.farg [6] * rayIn -> pcol [2]));
787 setcolor(talb, 0, 0, 0);
788 }
789 else {
790 copycolor(ctrans, cextinction);
791 copycolor(talb, salbedo);
792 }
793 }
794
795 else {
796 /* outside */
797 nratio = 1.0 / nratio;
798 setcolor(ctrans,
799 -mylog(mat -> oargs.farg [0] * rayIn -> pcol [0]),
800 -mylog(mat -> oargs.farg [1] * rayIn -> pcol [1]),
801 -mylog(mat -> oargs.farg [2] * rayIn -> pcol [2]));
802 setcolor(talb, 0, 0, 0);
803
804 if (mat -> otype == MAT_INTERFACE) {
805 setcolor(rayIn -> cext,
806 -mylog(mat -> oargs.farg [4] * rayIn -> pcol [0]),
807 -mylog(mat -> oargs.farg [5] * rayIn -> pcol [1]),
808 -mylog(mat -> oargs.farg [6] * rayIn -> pcol [2]));
809 setcolor(rayIn -> albedo, 0, 0, 0);
810 rayIn -> gecc = 0;
811 }
812 }
813
814 /* compute cos theta2 */
815 d2 = 1 - sqr(nratio) * (1 - sqr(cos1));
816
817 if (d2 < FTINY) {
818 /* Total reflection */
819 refl = cos2 = 1.0;
820 }
821 else {
822 /* Refraction, compute Fresnel's equations */
823 cos2 = sqrt(d2);
824 d1 = cos1;
825 d2 = nratio * cos2;
826 d1 = (d1 - d2) / (d1 + d2);
827 refl = sqr(d1);
828 d1 = 1 / cos1;
829 d2 = nratio / cos2;
830 d1 = (d1 - d2) / (d1 + d2);
831 refl += sqr(d1);
832 refl *= 0.5;
833 }
834
835 if (pmapRandom(rouletteState) > refl) {
836 /* Refraction */
837 photonRay(rayIn, &rayOut, PMAP_REFRACT, NULL);
838 d1 = nratio * cos1 - cos2;
839
840 for (i = 0; i < 3; i++)
841 rayOut.rdir [i] = nratio * rayIn -> rdir [i] + d1 * dnorm [i];
842
843 if (hastexture && DOT(rayOut.rdir, rayIn -> ron) * hastexture >= -FTINY) {
844 d1 *= hastexture;
845
846 for (i = 0; i < 3; i++)
847 rayOut.rdir [i] = nratio * rayIn -> rdir [i] +
848 d1 * rayIn -> ron [i];
849
850 normalize(rayOut.rdir);
851 }
852
853 copycolor(rayOut.cext, ctrans);
854 copycolor(rayOut.albedo, talb);
855 }
856
857 else {
858 /* Reflection */
859 photonRay(rayIn, &rayOut, PMAP_SPECREFL, NULL);
860 VSUM(rayOut.rdir, rayIn -> rdir, dnorm, 2 * cos1);
861
862 if (hastexture && DOT(rayOut.rdir, rayIn -> ron) * hastexture <= FTINY)
863 for (i = 0; i < 3; i++)
864 rayOut.rdir [i] = rayIn -> rdir [i] +
865 2 * rayIn -> rod * rayIn -> ron [i];
866 }
867
868 /* Ray is modified by medium defined by cext and albedo in
869 * photonParticipate() */
870 tracePhoton(&rayOut);
871
872 return 0;
873 }
874
875
876
877 static int glassPhotonScatter (OBJREC *mat, RAY *rayIn)
878 /* Generate new photon ray for glass material and recurse */
879 {
880 float albedo, xi, ptrans;
881 COLOR mcolor, refl, trans;
882 double pdot, cos2, d, r1e, r1m, rindex = 0.0;
883 FVECT pnorm, pdir;
884 int hastexture, i;
885 RAY rayOut;
886
887 /* check arguments */
888 if (mat -> oargs.nfargs == 3)
889 rindex = RINDEX;
890 else if (mat -> oargs.nfargs == 4)
891 rindex = mat -> oargs.farg [3];
892 else objerror(mat, USER, "bad arguments");
893
894 copycolor(mcolor, mat -> oargs.farg);
895
896 /* get modifiers */
897 raytexture(rayIn, mat -> omod);
898
899 /* reorient if necessary */
900 if (rayIn -> rod < 0)
901 flipsurface(rayIn);
902 if ((hastexture = DOT(rayIn -> pert, rayIn -> pert)) > FTINY * FTINY)
903 pdot = raynormal(pnorm, rayIn);
904 else {
905 VCOPY(pnorm, rayIn -> ron);
906 pdot = rayIn -> rod;
907 }
908
909 /* Modify material color */
910 multcolor(mcolor, rayIn -> pcol);
911
912 /* angular transmission */
913 cos2 = sqrt((1 - 1 / sqr(rindex)) + sqr(pdot / rindex));
914 setcolor(mcolor, pow(mcolor [0], 1 / cos2), pow(mcolor [1], 1 / cos2),
915 pow(mcolor [2], 1 / cos2));
916
917 /* compute reflection */
918 r1e = (pdot - rindex * cos2) / (pdot + rindex * cos2);
919 r1e *= r1e;
920 r1m = (1 / pdot - rindex / cos2) / (1 / pdot + rindex / cos2);
921 r1m *= r1m;
922
923 for (i = 0; i < 3; i++) {
924 double r1ed2, r1md2, d2;
925
926 d = mcolor [i];
927 d2 = sqr(d);
928 r1ed2 = sqr(r1e) * d2;
929 r1md2 = sqr(r1m) * d2;
930
931 /* compute transmittance */
932 trans [i] = 0.5 * d *
933 (sqr(1 - r1e) / (1 - r1ed2) + sqr(1 - r1m) / (1 - r1md2));
934
935 /* compute reflectance */
936 refl [i] = 0.5 * (r1e * (1 + (1 - 2 * r1e) * d2) / (1 - r1ed2) +
937 r1m * (1 + (1 - 2 * r1m) * d2) / (1 - r1md2));
938 }
939
940 /* Set up probabilities */
941 ptrans = colorAvg(trans);
942 albedo = colorAvg(refl) + ptrans;
943 xi = pmapRandom(rouletteState);
944
945
946 if (xi > albedo)
947 /* Absorbed */
948 return 0;
949
950 if (xi > (albedo -= ptrans)) {
951 /* Transmitted */
952
953 if (hastexture) {
954 /* perturb direction */
955 VSUM(pdir, rayIn -> rdir, rayIn -> pert, 2 * (1 - rindex));
956
957 if (normalize(pdir) == 0) {
958 objerror(mat, WARNING, "bad perturbation");
959 VCOPY(pdir, rayIn -> rdir);
960 }
961 }
962 else VCOPY(pdir, rayIn -> rdir);
963
964 VCOPY(rayOut.rdir, pdir);
965 photonRay(rayIn, &rayOut, PMAP_SPECTRANS, mcolor);
966 }
967
968 else {
969 /* reflected ray */
970 VSUM(rayOut.rdir, rayIn -> rdir, pnorm, 2 * pdot);
971 photonRay(rayIn, &rayOut, PMAP_SPECREFL, mcolor);
972 }
973
974 tracePhoton(&rayOut);
975 return 0;
976 }
977
978
979
980 static int aliasPhotonScatter (OBJREC *mat, RAY *rayIn)
981 /* Transfer photon scattering to alias target */
982 {
983 OBJECT aliasObj;
984 OBJREC aliasRec;
985
986 /* Straight replacement? */
987 if (!mat -> oargs.nsargs) {
988 mat = objptr(mat -> omod);
989 photonScatter [mat -> otype] (mat, rayIn);
990
991 return 0;
992 }
993
994 /* Else replace alias */
995 if (mat -> oargs.nsargs != 1)
996 objerror(mat, INTERNAL, "bad # string arguments");
997
998 aliasObj = lastmod(objndx(mat), mat -> oargs.sarg [0]);
999
1000 if (aliasObj < 0)
1001 objerror(mat, USER, "bad reference");
1002
1003 memcpy(&aliasRec, objptr(aliasObj), sizeof(OBJREC));
1004
1005 /* Substitute modifier */
1006 aliasRec.omod = mat -> omod;
1007
1008 /* Replacement scattering routine */
1009 photonScatter [aliasRec.otype] (&aliasRec, rayIn);
1010 return 0;
1011 }
1012
1013
1014
1015 static int clipPhotonScatter (OBJREC *mat, RAY *rayIn)
1016 /* Generate new photon ray for antimatter material and recurse */
1017 {
1018 OBJECT obj = objndx(mat), mod, cset [MAXSET + 1], *modset;
1019 int entering, inside = 0, i;
1020 const RAY *rp;
1021 RAY rayOut;
1022
1023 if ((modset = (OBJECT*)mat -> os) == NULL) {
1024 if (mat -> oargs.nsargs < 1 || mat -> oargs.nsargs > MAXSET)
1025 objerror(mat, USER, "bad # arguments");
1026
1027 modset = (OBJECT*)malloc((mat -> oargs.nsargs + 1) * sizeof(OBJECT));
1028
1029 if (modset == NULL)
1030 error(SYSTEM, "out of memory in clipPhotonScatter");
1031 modset [0] = 0;
1032
1033 for (i = 0; i < mat -> oargs.nsargs; i++) {
1034 if (!strcmp(mat -> oargs.sarg [i], VOIDID))
1035 continue;
1036
1037 if ((mod = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) {
1038 sprintf(errmsg, "unknown modifier \"%s\"", mat -> oargs.sarg [i]);
1039 objerror(mat, WARNING, errmsg);
1040 continue;
1041 }
1042
1043 if (inset(modset, mod)) {
1044 objerror(mat, WARNING, "duplicate modifier");
1045 continue;
1046 }
1047
1048 insertelem(modset, mod);
1049 }
1050
1051 mat -> os = (char*)modset;
1052 }
1053
1054 if (rayIn -> clipset != NULL)
1055 setcopy(cset, rayIn -> clipset);
1056 else cset [0] = 0;
1057
1058 entering = rayIn -> rod > 0;
1059
1060 /* Store photon incident from front if material defined as sensor */
1061 if (entering && inset(photonSensorSet, obj))
1062 addPhotons(rayIn);
1063
1064 for (i = modset [0]; i > 0; i--) {
1065 if (entering) {
1066 if (!inset(cset, modset [i])) {
1067 if (cset [0] >= MAXSET)
1068 error(INTERNAL, "set overflow in clipPhotonScatter");
1069 insertelem(cset, modset [i]);
1070 }
1071 }
1072 else if (inset(cset, modset [i]))
1073 deletelem(cset, modset [i]);
1074 }
1075
1076 rayIn -> newcset = cset;
1077
1078 if (strcmp(mat -> oargs.sarg [0], VOIDID)) {
1079 for (rp = rayIn; rp -> parent != NULL; rp = rp -> parent) {
1080 if ( !(rp -> rtype & RAYREFL) && rp->parent->ro != NULL &&
1081 inset(modset, rp -> parent -> ro -> omod)) {
1082
1083 if (rp -> parent -> rod > 0)
1084 inside++;
1085 else inside--;
1086 }
1087 }
1088
1089 if (inside > 0) {
1090 flipsurface(rayIn);
1091 mat = objptr(lastmod(obj, mat -> oargs.sarg [0]));
1092 photonScatter [mat -> otype] (mat, rayIn);
1093 return 0;
1094 }
1095 }
1096
1097 /* Else transfer ray */
1098 photonRay(rayIn, &rayOut, PMAP_XFER, NULL);
1099 tracePhoton(&rayOut);
1100
1101 return 0;
1102 }
1103
1104
1105
1106 static int mirrorPhotonScatter (OBJREC *mat, RAY *rayIn)
1107 /* Generate new photon ray for mirror material and recurse */
1108 {
1109 RAY rayOut;
1110 int rpure = 1, i;
1111 FVECT pnorm;
1112 double pdot;
1113 float albedo;
1114 COLOR mcolor;
1115
1116 /* check arguments */
1117 if (mat -> oargs.nfargs != 3 || mat -> oargs.nsargs > 1)
1118 objerror(mat, USER, "bad number of arguments");
1119
1120 /* back is black */
1121 if (rayIn -> rod < 0)
1122 return 0;
1123
1124 /* get modifiers */
1125 raytexture(rayIn, mat -> omod);
1126
1127 /* assign material color */
1128 copycolor(mcolor, mat -> oargs.farg);
1129 multcolor(mcolor, rayIn -> pcol);
1130
1131 /* Set up probabilities */
1132 albedo = colorAvg(mcolor);
1133
1134 if (pmapRandom(rouletteState) > albedo)
1135 /* Absorbed */
1136 return 0;
1137
1138 /* compute reflected ray */
1139 photonRay(rayIn, &rayOut, PMAP_SPECREFL, mcolor);
1140
1141 if (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)) {
1142 /* use textures */
1143 pdot = raynormal(pnorm, rayIn);
1144
1145 for (i = 0; i < 3; i++)
1146 rayOut.rdir [i] = rayIn -> rdir [i] + 2 * pdot * pnorm [i];
1147
1148 rpure = 0;
1149 }
1150
1151 /* Check for penetration */
1152 if (rpure || DOT(rayOut.rdir, rayIn -> ron) <= FTINY)
1153 for (i = 0; i < 3; i++)
1154 rayOut.rdir [i] = rayIn -> rdir [i] +
1155 2 * rayIn -> rod * rayIn -> ron [i];
1156
1157 tracePhoton(&rayOut);
1158 return 0;
1159 }
1160
1161
1162
1163 static int mistPhotonScatter (OBJREC *mat, RAY *rayIn)
1164 /* Generate new photon ray within mist and recurse */
1165 {
1166 COLOR mext;
1167 RREAL re, ge, be;
1168 RAY rayOut;
1169
1170 /* check arguments */
1171 if (mat -> oargs.nfargs > 7)
1172 objerror(mat, USER, "bad arguments");
1173
1174 if (mat -> oargs.nfargs > 2) {
1175 /* compute extinction */
1176 copycolor(mext, mat -> oargs.farg);
1177 /* get modifiers */
1178 raytexture(rayIn, mat -> omod);
1179 multcolor(mext, rayIn -> pcol);
1180 }
1181 else setcolor(mext, 0, 0, 0);
1182
1183 photonRay(rayIn, &rayOut, PMAP_XFER, NULL);
1184
1185 if (rayIn -> rod > 0) {
1186 /* entering ray */
1187 addcolor(rayOut.cext, mext);
1188
1189 if (mat -> oargs.nfargs > 5)
1190 copycolor(rayOut.albedo, mat -> oargs.farg + 3);
1191 if (mat -> oargs.nfargs > 6)
1192 rayOut.gecc = mat -> oargs.farg [6];
1193 }
1194
1195 else {
1196 /* leaving ray */
1197 re = max(rayIn -> cext [0] - mext [0], cextinction [0]);
1198 ge = max(rayIn -> cext [1] - mext [1], cextinction [1]);
1199 be = max(rayIn -> cext [2] - mext [2], cextinction [2]);
1200 setcolor(rayOut.cext, re, ge, be);
1201
1202 if (mat -> oargs.nfargs > 5)
1203 copycolor(rayOut.albedo, salbedo);
1204 if (mat -> oargs.nfargs > 6)
1205 rayOut.gecc = seccg;
1206 }
1207
1208 tracePhoton(&rayOut);
1209
1210 return 0;
1211 }
1212
1213
1214
1215 static int mx_dataPhotonScatter (OBJREC *mat, RAY *rayIn)
1216 /* Pass photon on to materials selected by mixture data */
1217 {
1218 OBJECT obj;
1219 double coef, pt [MAXDIM];
1220 DATARRAY *dp;
1221 OBJECT mod [2];
1222 MFUNC *mf;
1223 int i;
1224
1225 if (mat -> oargs.nsargs < 6)
1226 objerror(mat, USER, "bad # arguments");
1227
1228 obj = objndx(mat);
1229
1230 for (i = 0; i < 2; i++)
1231 if (!strcmp(mat -> oargs.sarg [i], VOIDID))
1232 mod [i] = OVOID;
1233 else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) {
1234 sprintf(errmsg, "undefined modifier \"%s\"", mat -> oargs.sarg [i]);
1235 objerror(mat, USER, errmsg);
1236 }
1237
1238 dp = getdata(mat -> oargs.sarg [3]);
1239 i = (1 << dp -> nd) - 1;
1240 mf = getfunc(mat, 4, i << 5, 0);
1241 setfunc(mat, rayIn);
1242 errno = 0;
1243
1244 for (i = 0; i < dp -> nd; i++) {
1245 pt [i] = evalue(mf -> ep [i]);
1246
1247 if (errno) {
1248 objerror(mat, WARNING, "compute error");
1249 return 0;
1250 }
1251 }
1252
1253 coef = datavalue(dp, pt);
1254 errno = 0;
1255 coef = funvalue(mat -> oargs.sarg [2], 1, &coef);
1256
1257 if (errno)
1258 objerror(mat, WARNING, "compute error");
1259 else {
1260 mat = objptr(mod [pmapRandom(rouletteState) < coef ? 0 : 1]);
1261 photonScatter [mat -> otype] (mat, rayIn);
1262 }
1263
1264 return 0;
1265 }
1266
1267
1268
1269 static int mx_pdataPhotonScatter (OBJREC *mat, RAY *rayIn)
1270 /* Pass photon on to materials selected by mixture picture */
1271 {
1272 OBJECT obj;
1273 double col [3], coef, pt [MAXDIM];
1274 DATARRAY *dp;
1275 OBJECT mod [2];
1276 MFUNC *mf;
1277 int i;
1278
1279 if (mat -> oargs.nsargs < 7)
1280 objerror(mat, USER, "bad # arguments");
1281
1282 obj = objndx(mat);
1283
1284 for (i = 0; i < 2; i++)
1285 if (!strcmp(mat -> oargs.sarg [i], VOIDID))
1286 mod [i] = OVOID;
1287 else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) {
1288 sprintf(errmsg, "undefined modifier \"%s\"", mat -> oargs.sarg [i]);
1289 objerror(mat, USER, errmsg);
1290 }
1291
1292 dp = getpict(mat -> oargs.sarg [3]);
1293 mf = getfunc(mat, 4, 0x3 << 5, 0);
1294 setfunc(mat, rayIn);
1295 errno = 0;
1296 pt [1] = evalue(mf -> ep [0]);
1297 pt [0] = evalue(mf -> ep [1]);
1298
1299 if (errno) {
1300 objerror(mat, WARNING, "compute error");
1301 return 0;
1302 }
1303
1304 for (i = 0; i < 3; i++)
1305 col [i] = datavalue(dp + i, pt);
1306
1307 errno = 0;
1308 coef = funvalue(mat -> oargs.sarg [2], 3, col);
1309
1310 if (errno)
1311 objerror(mat, WARNING, "compute error");
1312 else {
1313 mat = objptr(mod [pmapRandom(rouletteState) < coef ? 0 : 1]);
1314 photonScatter [mat -> otype] (mat, rayIn);
1315 }
1316
1317 return 0;
1318 }
1319
1320
1321
1322 static int mx_funcPhotonScatter (OBJREC *mat, RAY *rayIn)
1323 /* Pass photon on to materials selected by mixture function */
1324 {
1325 OBJECT obj, mod [2];
1326 int i;
1327 double coef;
1328 MFUNC *mf;
1329
1330 if (mat -> oargs.nsargs < 4)
1331 objerror(mat, USER, "bad # arguments");
1332
1333 obj = objndx(mat);
1334
1335 for (i = 0; i < 2; i++)
1336 if (!strcmp(mat -> oargs.sarg [i], VOIDID))
1337 mod [i] = OVOID;
1338 else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) {
1339 sprintf(errmsg, "undefined modifier \"%s\"", mat -> oargs.sarg [i]);
1340 objerror(mat, USER, errmsg);
1341 }
1342
1343 mf = getfunc(mat, 3, 0x4, 0);
1344 setfunc(mat, rayIn);
1345 errno = 0;
1346
1347 /* bound coefficient */
1348 coef = min(1, max(0, evalue(mf -> ep [0])));
1349
1350 if (errno)
1351 objerror(mat, WARNING, "compute error");
1352 else {
1353 mat = objptr(mod [pmapRandom(rouletteState) < coef ? 0 : 1]);
1354 photonScatter [mat -> otype] (mat, rayIn);
1355 }
1356
1357 return 0;
1358 }
1359
1360
1361
1362 static int pattexPhotonScatter (OBJREC *mat, RAY *rayIn)
1363 /* Generate new photon ray for pattern or texture modifier and recurse.
1364 This code is brought to you by Henkel! :^) */
1365 {
1366 RAY rayOut;
1367
1368 /* Get pattern */
1369 ofun [mat -> otype].funp(mat, rayIn);
1370 if (mat -> omod != OVOID) {
1371 /* Scatter using modifier (if any) */
1372 mat = objptr(mat -> omod);
1373 photonScatter [mat -> otype] (mat, rayIn);
1374 }
1375 else {
1376 /* Transfer ray if no modifier */
1377 photonRay(rayIn, &rayOut, PMAP_XFER, NULL);
1378 tracePhoton(&rayOut);
1379 }
1380
1381 return 0;
1382 }
1383
1384
1385
1386 #if 0
1387 static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn)
1388 /* Generate new photon ray for BSDF modifier and recurse. */
1389 {
1390 int hitFront;
1391 SDError err;
1392 FVECT upvec;
1393 MFUNC *mf;
1394 BSDFDAT nd;
1395 RAY rayOut;
1396
1397 /* Following code adapted from m_bsdf() */
1398 /* Check arguments */
1399 if (mat -> oargs.nsargs < 6 || mat -> oargs.nfargs > 9 ||
1400 mat -> oargs.nfargs % 3)
1401 objerror(mat, USER, "bad # arguments");
1402
1403 hitFront = (rayIn -> rod > 0);
1404
1405 /* Load cal file */
1406 mf = getfunc(mat, 5, 0x1d, 1);
1407
1408 /* Get thickness */
1409 nd.thick = evalue(mf -> ep [0]);
1410 if ((-FTINY <= nd.thick) & (nd.thick <= FTINY))
1411 nd.thick = .0;
1412
1413 if (nd.thick != .0 || (!hitFront && !backvis)) {
1414 /* Proxy geometry present, so use it instead and transfer ray */
1415 photonRay(rayIn, &rayOut, PMAP_XFER, NULL);
1416 tracePhoton(&rayOut);
1417
1418 return 0;
1419 }
1420
1421 /* Get BSDF data */
1422 nd.sd = loadBSDF(mat -> oargs.sarg [1]);
1423
1424 /* Diffuse reflectance */
1425 if (hitFront) {
1426 if (mat -> oargs.nfargs < 3)
1427 setcolor(nd.rdiff, .0, .0, .0);
1428 else setcolor(nd.rdiff, mat -> oargs.farg [0], mat -> oargs.farg [1],
1429 mat -> oargs.farg [2]);
1430 }
1431 else if (mat -> oargs.nfargs < 6) {
1432 /* Check for absorbing backside */
1433 if (!backvis && !nd.sd -> rb && !nd.sd -> tf) {
1434 SDfreeCache(nd.sd);
1435 return 0;
1436 }
1437
1438 setcolor(nd.rdiff, .0, .0, .0);
1439 }
1440 else setcolor(nd.rdiff, mat -> oargs.farg [3], mat -> oargs.farg [4],
1441 mat -> oargs.farg [5]);
1442
1443 /* Diffuse transmittance */
1444 if (mat -> oargs.nfargs < 9)
1445 setcolor(nd.tdiff, .0, .0, .0);
1446 else setcolor(nd.tdiff, mat -> oargs.farg [6], mat -> oargs.farg [7],
1447 mat -> oargs.farg [8]);
1448
1449 nd.mp = mat;
1450 nd.pr = rayIn;
1451
1452 /* Get modifiers */
1453 raytexture(rayIn, mat -> omod);
1454
1455 /* Modify diffuse values */
1456 multcolor(nd.rdiff, rayIn -> pcol);
1457 multcolor(nd.tdiff, rayIn -> pcol);
1458
1459 /* Get up vector & xform to world coords */
1460 upvec [0] = evalue(mf -> ep [1]);
1461 upvec [1] = evalue(mf -> ep [2]);
1462 upvec [2] = evalue(mf -> ep [3]);
1463
1464 if (mf -> fxp != &unitxf) {
1465 multv3(upvec, upvec, mf -> fxp -> xfm);
1466 nd.thick *= mf -> fxp -> sca;
1467 }
1468
1469 if (rayIn -> rox) {
1470 multv3(upvec, upvec, rayIn -> rox -> f.xfm);
1471 nd.thick *= rayIn -> rox -> f.sca;
1472 }
1473
1474 /* Perturb normal */
1475 raynormal(nd.pnorm, rayIn);
1476
1477 /* Xform incident dir to local BSDF coords */
1478 err = SDcompXform(nd.toloc, nd.pnorm, upvec);
1479
1480 if (!err) {
1481 nd.vray [0] = -rayIn -> rdir [0];
1482 nd.vray [1] = -rayIn -> rdir [1];
1483 nd.vray [2] = -rayIn -> rdir [2];
1484 err = SDmapDir(nd.vray, nd.toloc, nd.vray);
1485 }
1486
1487 if (!err)
1488 err = SDinvXform(nd.fromloc, nd.toloc);
1489
1490 if (err) {
1491 objerror(mat, WARNING, "Illegal orientation vector");
1492 return 0;
1493 }
1494
1495 /* Determine BSDF resolution */
1496 err = SDsizeBSDF(nd.sr_vpsa, nd.vray, NULL, SDqueryMin + SDqueryMax, nd.sd);
1497
1498 if (err)
1499 objerror(mat, USER, transSDError(err));
1500
1501 nd.sr_vpsa [0] = sqrt(nd.sr_vpsa [0]);
1502 nd.sr_vpsa [1] = sqrt(nd.sr_vpsa [1]);
1503
1504 /* Orient perturbed normal towards incident side */
1505 if (!hitFront) {
1506 nd.pnorm [0] = -nd.pnorm [0];
1507 nd.pnorm [1] = -nd.pnorm [1];
1508 nd.pnorm [2] = -nd.pnorm [2];
1509 }
1510
1511 /* Following code adapted from SDsampBSDF() */
1512 {
1513 SDSpectralDF *rdf, *tdf;
1514 SDValue bsdfVal;
1515 double xi, rhoDiff = 0;
1516 float coef [SDmaxCh];
1517 int i, j, n, nr;
1518 SDComponent *sdc;
1519 const SDCDst **cdarr = NULL;
1520
1521 /* Get diffuse albedo (?) */
1522 if (hitFront) {
1523 bsdfVal = nd.sd -> rLambFront;
1524 rdf = nd.sd -> rf;
1525 tdf = nd.sd -> tf ? nd.sd -> tf : nd.sd -> tb;
1526 }
1527 else {
1528 bsdfVal = nd.sd -> rLambBack;
1529 rdf = nd.sd -> rb;
1530 tdf = nd.sd -> tb ? nd.sd -> tb : nd.sd -> tf;
1531 }
1532
1533 rhoDiff = bsdfVal.cieY;
1534 bsdfVal.cieY += nd.sd -> tLamb.cieY;
1535
1536 /* Allocate non-diffuse sampling */
1537 i = nr = rdf ? rdf -> ncomp : 0;
1538 j = tdf ? tdf -> ncomp : 0;
1539 n = i + j;
1540
1541 if (n > 0 && !(cdarr = (const SDCDst**)malloc(n * sizeof(SDCDst*))))
1542 objerror(mat, USER, transSDError(SDEmemory));
1543
1544 while (j-- > 0) {
1545 /* Sum up non-diffuse transmittance */
1546 cdarr [i + j] = (*tdf -> comp [j].func -> getCDist)(nd.vray, &tdf -> comp [j]);
1547
1548 if (!cdarr [i + j])
1549 cdarr [i + j] = &SDemptyCD;
1550 else bsdfVal.cieY += cdarr [i + j] -> cTotal;
1551 }
1552
1553 while (i-- > 0) {
1554 /* Sum up non-diffuse reflectance */
1555 cdarr [i] = (*rdf -> comp [i].func -> getCDist)(nd.vray, &rdf -> comp [i]);
1556
1557 if (!cdarr [i])
1558 cdarr [i] = &SDemptyCD;
1559 else bsdfVal.cieY += cdarr [i] -> cTotal;
1560 }
1561
1562 if (bsdfVal.cieY <= FTINY) {
1563 /* Don't bother sampling, just absorb photon */
1564 if (cdarr)
1565 free(cdarr);
1566 return 0;
1567 }
1568
1569 /* Insert direct and indirect photon hits if diffuse component */
1570 if (rhoDiff > FTINY || nd.sd -> tLamb.cieY > FTINY)
1571 addPhotons(rayIn);
1572
1573 xi = pmapRandom(rouletteState);
1574
1575 if ((xi -= rhoDiff) <= 0) {
1576 /* Diffuse reflection */
1577 photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff);
1578 diffPhotonScatter(nd.pnorm, &rayOut);
1579 }
1580 else if ((xi -= nd.sd -> tLamb.cieY) <= 0) {
1581 /* Diffuse transmission */
1582 flipsurface(rayIn);
1583 photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff);
1584 bsdfVal.spec = nd.sd -> tLamb.spec;
1585 diffPhotonScatter(nd.pnorm, &rayOut);
1586 }
1587 else {
1588 int rayOutType;
1589 COLOR bsdfRGB;
1590
1591 /* Non-diffuse CDF inversion (?) */
1592 for (i = 0; i < n && (xi -= cdarr [i] -> cTotal) > 0; i++);
1593
1594 if (i >= n) {
1595 /* Absorbed -- photon went Deer Hunter */
1596 if (cdarr)
1597 free(cdarr);
1598 return 0;
1599 }
1600
1601 if (i < nr) {
1602 /* Non-diffuse reflection */
1603 sdc = &rdf -> comp [i];
1604 rayOutType = PMAP_SPECREFL;
1605 }
1606 else {
1607 /* Non-diffuse transmission */
1608 sdc = &tdf -> comp [i - nr];
1609 rayOutType = PMAP_SPECTRANS;
1610 }
1611
1612 /* Generate non-diff sample dir */
1613 VCOPY(rayOut.rdir, nd.vray);
1614 err = (*sdc -> func -> sampCDist)
1615 (rayOut.rdir, pmapRandom(scatterState), cdarr [i]);
1616 if (err)
1617 objerror(mat, USER, transSDError(SDEinternal));
1618
1619 /* Get colour */
1620 j = (*sdc -> func -> getBSDFs)(coef, rayOut.rdir, nd.vray, sdc);
1621
1622 if (j <= 0) {
1623 sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error",
1624 nd.sd -> name);
1625 objerror(mat, USER, transSDError(SDEinternal));
1626 }
1627
1628 bsdfVal.spec = sdc -> cspec [0];
1629 rhoDiff = coef [0];
1630
1631 while (--j) {
1632 c_cmix(&bsdfVal.spec, rhoDiff, &bsdfVal.spec, coef [j],
1633 &sdc -> cspec [j]);
1634 rhoDiff += coef [j];
1635 }
1636
1637 /* ? */
1638 c_ccvt(&bsdfVal.spec, C_CSXY + C_CSSPEC);
1639 ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB);
1640
1641 /* Xform outgoing dir to world coords */
1642 if ((err = SDmapDir(rayOut.rdir, nd.fromloc, rayOut.rdir))) {
1643 objerror(mat, USER, transSDError(err));
1644 return 0;
1645 }
1646
1647 photonRay(rayIn, &rayOut, rayOutType, bsdfRGB);
1648 }
1649
1650 if (cdarr)
1651 free(cdarr);
1652 }
1653
1654 /* Clean up BSDF */
1655 SDfreeCache(nd.sd);
1656
1657 tracePhoton(&rayOut);
1658 return 0;
1659 }
1660 #else
1661
1662 static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn)
1663 /* Generate new photon ray for BSDF modifier and recurse. */
1664 {
1665 int hitFront;
1666 SDError err;
1667 SDValue bsdfVal;
1668 FVECT upvec;
1669 MFUNC *mf;
1670 BSDFDAT nd;
1671 RAY rayOut;
1672 COLOR bsdfRGB;
1673 double prDiff, ptDiff, prDiffSD, ptDiffSD, prSpecSD, ptSpecSD,
1674 albedo, xi, xi2;
1675 const double patAlb = colorAvg(rayIn -> pcol);
1676
1677 /* Following code adapted from m_bsdf() */
1678 /* Check arguments */
1679 if (mat -> oargs.nsargs < 6 || mat -> oargs.nfargs > 9 ||
1680 mat -> oargs.nfargs % 3)
1681 objerror(mat, USER, "bad # arguments");
1682
1683 hitFront = (rayIn -> rod > 0);
1684
1685 /* Load cal file */
1686 mf = getfunc(mat, 5, 0x1d, 1);
1687
1688 /* Get thickness */
1689 nd.thick = evalue(mf -> ep [0]);
1690 if ((-FTINY <= nd.thick) & (nd.thick <= FTINY))
1691 nd.thick = .0;
1692
1693 if (nd.thick != .0 || (!hitFront && !backvis)) {
1694 /* Proxy geometry present, so use it instead and transfer ray */
1695 photonRay(rayIn, &rayOut, PMAP_XFER, NULL);
1696 tracePhoton(&rayOut);
1697
1698 return 0;
1699 }
1700
1701 /* Get BSDF data */
1702 nd.sd = loadBSDF(mat -> oargs.sarg [1]);
1703
1704 /* Extra diffuse reflectance from material def */
1705 if (hitFront) {
1706 if (mat -> oargs.nfargs < 3)
1707 setcolor(nd.rdiff, .0, .0, .0);
1708 else setcolor(nd.rdiff, mat -> oargs.farg [0], mat -> oargs.farg [1],
1709 mat -> oargs.farg [2]);
1710 }
1711 else if (mat -> oargs.nfargs < 6) {
1712 /* Check for absorbing backside */
1713 if (!backvis && !nd.sd -> rb && !nd.sd -> tf) {
1714 SDfreeCache(nd.sd);
1715 return 0;
1716 }
1717
1718 setcolor(nd.rdiff, .0, .0, .0);
1719 }
1720 else setcolor(nd.rdiff, mat -> oargs.farg [3], mat -> oargs.farg [4],
1721 mat -> oargs.farg [5]);
1722
1723 /* Extra diffuse transmittance from material def */
1724 if (mat -> oargs.nfargs < 9)
1725 setcolor(nd.tdiff, .0, .0, .0);
1726 else setcolor(nd.tdiff, mat -> oargs.farg [6], mat -> oargs.farg [7],
1727 mat -> oargs.farg [8]);
1728
1729 nd.mp = mat;
1730 nd.pr = rayIn;
1731
1732 /* Get modifiers */
1733 raytexture(rayIn, mat -> omod);
1734
1735 /* Modify diffuse values */
1736 multcolor(nd.rdiff, rayIn -> pcol);
1737 multcolor(nd.tdiff, rayIn -> pcol);
1738
1739 /* Get up vector & xform to world coords */
1740 upvec [0] = evalue(mf -> ep [1]);
1741 upvec [1] = evalue(mf -> ep [2]);
1742 upvec [2] = evalue(mf -> ep [3]);
1743
1744 if (mf -> fxp != &unitxf) {
1745 multv3(upvec, upvec, mf -> fxp -> xfm);
1746 nd.thick *= mf -> fxp -> sca;
1747 }
1748
1749 if (rayIn -> rox) {
1750 multv3(upvec, upvec, rayIn -> rox -> f.xfm);
1751 nd.thick *= rayIn -> rox -> f.sca;
1752 }
1753
1754 /* Perturb normal */
1755 raynormal(nd.pnorm, rayIn);
1756
1757 /* Xform incident dir to local BSDF coords */
1758 err = SDcompXform(nd.toloc, nd.pnorm, upvec);
1759
1760 if (!err) {
1761 nd.vray [0] = -rayIn -> rdir [0];
1762 nd.vray [1] = -rayIn -> rdir [1];
1763 nd.vray [2] = -rayIn -> rdir [2];
1764 err = SDmapDir(nd.vray, nd.toloc, nd.vray);
1765 }
1766
1767 if (!err)
1768 err = SDinvXform(nd.fromloc, nd.toloc);
1769
1770 if (err) {
1771 objerror(mat, WARNING, "Illegal orientation vector");
1772 return 0;
1773 }
1774
1775 /* Determine BSDF resolution */
1776 err = SDsizeBSDF(nd.sr_vpsa, nd.vray, NULL, SDqueryMin + SDqueryMax, nd.sd);
1777
1778 if (err)
1779 objerror(mat, USER, transSDError(err));
1780
1781 nd.sr_vpsa [0] = sqrt(nd.sr_vpsa [0]);
1782 nd.sr_vpsa [1] = sqrt(nd.sr_vpsa [1]);
1783
1784 /* Orient perturbed normal towards incident side */
1785 if (!hitFront) {
1786 nd.pnorm [0] = -nd.pnorm [0];
1787 nd.pnorm [1] = -nd.pnorm [1];
1788 nd.pnorm [2] = -nd.pnorm [2];
1789 }
1790
1791 /* Get scatter probabilities (weighted by pattern except for spec refl)
1792 * prDiff, ptDiff: extra diffuse component in material def
1793 * prDiffSD, ptDiffSD: diffuse (constant) component in SDF
1794 * prSpecSD, ptSpecSD: non-diffuse ("specular") component in SDF
1795 * albedo: sum of above, inverse absorption probability */
1796 prDiff = colorAvg(nd.rdiff);
1797 ptDiff = colorAvg(nd.tdiff);
1798 prDiffSD = patAlb * SDdirectHemi(nd.vray, SDsampDf | SDsampR, nd.sd);
1799 ptDiffSD = patAlb * SDdirectHemi(nd.vray, SDsampDf | SDsampT, nd.sd);
1800 prSpecSD = SDdirectHemi(nd.vray, SDsampSp | SDsampR, nd.sd);
1801 ptSpecSD = patAlb * SDdirectHemi(nd.vray, SDsampSp | SDsampT, nd.sd);
1802 albedo = prDiff + ptDiff + prDiffSD + ptDiffSD + prSpecSD + ptSpecSD;
1803
1804 /*
1805 if (albedo > 1)
1806 objerror(mat, WARNING, "Invalid albedo");
1807 */
1808
1809 /* Insert direct and indirect photon hits if diffuse component */
1810 if (prDiff + ptDiff + prDiffSD + ptDiffSD > FTINY)
1811 addPhotons(rayIn);
1812
1813 xi = xi2 = pmapRandom(rouletteState);
1814
1815 if (xi > albedo)
1816 /* Absorbtion */
1817 return 0;
1818
1819 if ((xi -= prDiff) <= 0) {
1820 /* Diffuse reflection (extra component in material def) */
1821 photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff);
1822 diffPhotonScatter(nd.pnorm, &rayOut);
1823 }
1824
1825 else if ((xi -= ptDiff) <= 0) {
1826 /* Diffuse transmission (extra component in material def) */
1827 flipsurface(rayIn);
1828 photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff);
1829 diffPhotonScatter(nd.pnorm, &rayOut);
1830 }
1831
1832 else { /* Sample SDF */
1833 if ((xi -= prDiffSD) <= 0) {
1834 /* Diffuse SDF reflection (constant component) */
1835 if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2,
1836 SDsampDf | SDsampR, nd.sd)))
1837 objerror(mat, USER, transSDError(err));
1838
1839 /* Apply pattern to spectral component */
1840 ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB);
1841 multcolor(bsdfRGB, rayIn -> pcol);
1842 photonRay(rayIn, &rayOut, PMAP_DIFFREFL, bsdfRGB);
1843 }
1844
1845 else if ((xi -= ptDiffSD) <= 0) {
1846 /* Diffuse SDF transmission (constant component) */
1847 if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2,
1848 SDsampDf | SDsampT, nd.sd)))
1849 objerror(mat, USER, transSDError(err));
1850
1851 /* Apply pattern to spectral component */
1852 ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB);
1853 multcolor(bsdfRGB, rayIn -> pcol);
1854 addcolor(bsdfRGB, nd.tdiff);
1855 flipsurface(rayIn); /* Necessary? */
1856 photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, bsdfRGB);
1857 }
1858
1859 else if ((xi -= prSpecSD) <= 0) {
1860 /* Non-diffuse ("specular") SDF reflection */
1861 if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2,
1862 SDsampSp | SDsampR, nd.sd)))
1863 objerror(mat, USER, transSDError(err));
1864
1865 ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB);
1866 photonRay(rayIn, &rayOut, PMAP_SPECREFL, bsdfRGB);
1867 }
1868
1869 else {
1870 /* Non-diffuse ("specular") SDF transmission */
1871 if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2,
1872 SDsampSp | SDsampT, nd.sd)))
1873 objerror(mat, USER, transSDError(err));
1874
1875 /* Apply pattern to spectral component */
1876 ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB);
1877 multcolor(bsdfRGB, rayIn -> pcol);
1878 flipsurface(rayIn); /* Necessary? */
1879 photonRay(rayIn, &rayOut, PMAP_SPECTRANS, bsdfRGB);
1880 }
1881
1882 /* Xform outgoing dir to world coords */
1883 if ((err = SDmapDir(rayOut.rdir, nd.fromloc, nd.vray))) {
1884 objerror(mat, USER, transSDError(err));
1885 return 0;
1886 }
1887 }
1888
1889 /* Clean up */
1890 SDfreeCache(nd.sd);
1891
1892 tracePhoton(&rayOut);
1893 return 0;
1894 }
1895 #endif
1896
1897
1898
1899 static int lightPhotonScatter (OBJREC* mat, RAY* ray)
1900 /* Light sources doan' reflect */
1901 {
1902 return 0;
1903 }
1904
1905
1906
1907 void initPhotonScatterFuncs ()
1908 /* Init photonScatter[] dispatch table */
1909 {
1910 int i;
1911
1912 for (i = 0; i < NUMOTYPE; i++)
1913 photonScatter [i] = o_default;
1914
1915 photonScatter [MAT_LIGHT] = photonScatter [MAT_ILLUM] =
1916 photonScatter [MAT_GLOW] = photonScatter [MAT_SPOT] =
1917 lightPhotonScatter;
1918
1919 photonScatter [MAT_PLASTIC] = photonScatter [MAT_METAL] =
1920 photonScatter [MAT_TRANS] = normalPhotonScatter;
1921
1922 photonScatter [MAT_PLASTIC2] = photonScatter [MAT_METAL2] =
1923 photonScatter [MAT_TRANS2] = anisoPhotonScatter;
1924
1925 photonScatter [MAT_DIELECTRIC] = photonScatter [MAT_INTERFACE] =
1926 dielectricPhotonScatter;
1927
1928 photonScatter [MAT_MIST] = mistPhotonScatter;
1929 photonScatter [MAT_GLASS] = glassPhotonScatter;
1930 photonScatter [MAT_CLIP] = clipPhotonScatter;
1931 photonScatter [MAT_MIRROR] = mirrorPhotonScatter;
1932 photonScatter [MIX_FUNC] = mx_funcPhotonScatter;
1933 photonScatter [MIX_DATA] = mx_dataPhotonScatter;
1934 photonScatter [MIX_PICT]= mx_pdataPhotonScatter;
1935
1936 photonScatter [PAT_BDATA] = photonScatter [PAT_CDATA] =
1937 photonScatter [PAT_BFUNC] = photonScatter [PAT_CFUNC] =
1938 photonScatter [PAT_CPICT] = photonScatter [TEX_FUNC] =
1939 photonScatter [TEX_DATA] = pattexPhotonScatter;
1940
1941 photonScatter [MOD_ALIAS] = aliasPhotonScatter;
1942 photonScatter [MAT_BSDF] = bsdfPhotonScatter;
1943 }