ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapmat.c
Revision: 2.5
Committed: Thu May 21 05:54:54 2015 UTC (9 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +3 -3 lines
Log Message:
Made axis randomization optional in getperpendicular()

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