ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapmat.c
Revision: 2.17
Committed: Mon Jun 25 20:49:10 2018 UTC (5 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.16: +16 -11 lines
Log Message:
Added sBSDF material type for explicit control over peak extraction

File Contents

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