ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapmat.c
Revision: 2.3
Committed: Fri May 8 13:20:23 2015 UTC (9 years ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.2: +9 -2 lines
Log Message:
Double-counting bugfix for glow sources (thanks DGM!), revised copyright

File Contents

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