ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapmat.c
Revision: 2.8
Committed: Tue Aug 18 18:45:55 2015 UTC (8 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +4 -1 lines
Log Message:
Added missing RCSid forgotten during initial check-in

File Contents

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