ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapmat.c
Revision: 2.12
Committed: Tue Feb 23 12:42:41 2016 UTC (8 years, 3 months ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.11: +1 -1 lines
Log Message:
Revised PMAP overcounting macros in pmapmat.h, added one in badcomponent()
for overcounted caustics from transferred shadow/ambient rays

File Contents

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