ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapsrc.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: +3 -2 lines
Log Message:
Double-counting bugfix for glow sources (thanks DGM!), revised copyright

File Contents

# User Rev Content
1 greg 2.1 /*
2     ==================================================================
3     Photon map support routines for emission from light sources
4    
5     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
6     (c) Fraunhofer Institute for Solar Energy Systems,
7 rschregle 2.3 (c) Lucerne University of Applied Sciences and Arts,
8     supported by the Swiss National Science Foundation (SNSF, #147053)
9 greg 2.1 ==================================================================
10    
11 rschregle 2.3 $Id: pmapsrc.c,v 2.2 2015/04/21 19:16:51 greg Exp $
12 greg 2.1 */
13    
14    
15    
16     #include "pmapsrc.h"
17     #include "pmap.h"
18     #include "pmaprand.h"
19     #include "otypes.h"
20    
21    
22    
23     SRCREC *photonPorts = NULL; /* Photon port list */
24     unsigned numPhotonPorts = 0;
25    
26     void (*photonPartition [NUMOTYPE]) (EmissionMap*);
27     void (*photonOrigin [NUMOTYPE]) (EmissionMap*);
28    
29     extern OBJECT ambset [];
30    
31    
32    
33     static int flatPhotonPartition2 (EmissionMap* emap, unsigned long mp,
34     FVECT cent, FVECT u, FVECT v,
35     double du2, double dv2)
36     /* Recursive part of flatPhotonPartition(..) */
37     {
38     FVECT newct, newax;
39     unsigned long npl, npu;
40    
41     if (mp > emap -> maxPartitions) {
42     /* Enlarge partition array */
43     emap -> maxPartitions <<= 1;
44     emap -> partitions = (unsigned char*)realloc(emap -> partitions,
45     emap -> maxPartitions >> 1);
46    
47     if (!emap -> partitions)
48     error(USER, "can't allocate source partitions");
49    
50 greg 2.2 memset(emap -> partitions + (emap -> maxPartitions >> 2), 0,
51     emap -> maxPartitions >> 2);
52 greg 2.1 }
53    
54     if (du2 * dv2 <= 1) { /* hit limit? */
55     setpart(emap -> partitions, emap -> partitionCnt, S0);
56     emap -> partitionCnt++;
57     return 1;
58     }
59    
60     if (du2 > dv2) { /* subdivide in U */
61     setpart(emap -> partitions, emap -> partitionCnt, SU);
62     emap -> partitionCnt++;
63     newax [0] = 0.5 * u [0];
64     newax [1] = 0.5 * u [1];
65     newax [2] = 0.5 * u [2];
66     u = newax;
67     du2 *= 0.25;
68     }
69    
70     else { /* subdivide in V */
71     setpart(emap -> partitions, emap -> partitionCnt, SV);
72     emap -> partitionCnt++;
73     newax [0] = 0.5 * v [0];
74     newax [1] = 0.5 * v [1];
75     newax [2] = 0.5 * v [2];
76     v = newax;
77     dv2 *= 0.25;
78     }
79    
80     /* lower half */
81     newct [0] = cent [0] - newax [0];
82     newct [1] = cent [1] - newax [1];
83     newct [2] = cent [2] - newax [2];
84     npl = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2);
85    
86     /* upper half */
87     newct [0] = cent [0] + newax [0];
88     newct [1] = cent [1] + newax [1];
89     newct [2] = cent [2] + newax [2];
90     npu = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2);
91    
92     /* return total */
93     return npl + npu;
94     }
95    
96    
97    
98     static void flatPhotonPartition (EmissionMap* emap)
99     /* Partition flat source for photon emission */
100     {
101     RREAL *vp;
102     double du2, dv2;
103    
104 greg 2.2 memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
105 greg 2.1 emap -> partArea = srcsizerat * thescene.cusize;
106     emap -> partArea *= emap -> partArea;
107     vp = emap -> src -> ss [SU];
108     du2 = DOT(vp, vp) / emap -> partArea;
109     vp = emap -> src -> ss [SV];
110     dv2 = DOT(vp, vp) / emap -> partArea;
111     emap -> partitionCnt = 0;
112     emap -> numPartitions = flatPhotonPartition2(emap, 1, emap -> src -> sloc,
113     emap -> src -> ss [SU],
114     emap -> src -> ss [SV],
115     du2, dv2);
116     emap -> partitionCnt = 0;
117     emap -> partArea = emap -> src -> ss2 / emap -> numPartitions;
118     }
119    
120    
121    
122     static void sourcePhotonPartition (EmissionMap* emap)
123     /* Partition scene cube faces or photon port for photon emission from
124     distant source */
125     {
126     if (emap -> port) {
127     /* Partition photon port */
128     SRCREC *src = emap -> src;
129     emap -> src = emap -> port;
130     photonPartition [emap -> src -> so -> otype] (emap);
131     emap -> src = src;
132     }
133    
134     else {
135     /* No photon ports defined, so partition scene cube faces */
136 greg 2.2 memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
137 greg 2.1 setpart(emap -> partitions, 0, S0);
138     emap -> partitionCnt = 0;
139     emap -> numPartitions = 1 / srcsizerat;
140     emap -> numPartitions *= emap -> numPartitions;
141     emap -> partArea = sqr(thescene.cusize) / emap -> numPartitions;
142     emap -> numPartitions *= 6;
143     }
144     }
145    
146    
147    
148     static void spherePhotonPartition (EmissionMap* emap)
149     /* Partition spherical source into equal solid angles using uniform
150     mapping for photon emission */
151     {
152     unsigned numTheta, numPhi;
153    
154 greg 2.2 memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
155 greg 2.1 setpart(emap -> partitions, 0, S0);
156     emap -> partArea = 4 * PI * sqr(emap -> src -> srad);
157     emap -> numPartitions = emap -> partArea /
158     sqr(srcsizerat * thescene.cusize);
159    
160     numTheta = max(sqrt(2 * emap -> numPartitions / PI) + 0.5, 1);
161     numPhi = 0.5 * PI * numTheta + 0.5;
162    
163     emap -> numPartitions = (unsigned long)numTheta * numPhi;
164     emap -> partitionCnt = 0;
165     emap -> partArea /= emap -> numPartitions;
166     }
167    
168    
169    
170     static int cylPhotonPartition2 (EmissionMap* emap, unsigned long mp,
171     FVECT cent, FVECT axis, double d2)
172     /* Recursive part of cyPhotonPartition(..) */
173     {
174     FVECT newct, newax;
175     unsigned long npl, npu;
176    
177     if (mp > emap -> maxPartitions) {
178     /* Enlarge partition array */
179     emap -> maxPartitions <<= 1;
180     emap -> partitions = (unsigned char*)realloc(emap -> partitions,
181     emap -> maxPartitions >> 1);
182     if (!emap -> partitions)
183     error(USER, "can't allocate source partitions");
184 greg 2.2
185     memset(emap -> partitions + (emap -> maxPartitions >> 2), 0,
186 greg 2.1 emap -> maxPartitions >> 2);
187     }
188    
189     if (d2 <= 1) {
190     /* hit limit? */
191     setpart(emap -> partitions, emap -> partitionCnt, S0);
192     emap -> partitionCnt++;
193     return 1;
194     }
195    
196     /* subdivide */
197     setpart(emap -> partitions, emap -> partitionCnt, SU);
198     emap -> partitionCnt++;
199     newax [0] = 0.5 * axis [0];
200     newax [1] = 0.5 * axis [1];
201     newax [2] = 0.5 * axis [2];
202     d2 *= 0.25;
203    
204     /* lower half */
205     newct [0] = cent [0] - newax [0];
206     newct [1] = cent [1] - newax [1];
207     newct [2] = cent [2] - newax [2];
208     npl = cylPhotonPartition2(emap, mp << 1, newct, newax, d2);
209    
210     /* upper half */
211     newct [0] = cent [0] + newax [0];
212     newct [1] = cent [1] + newax [1];
213     newct [2] = cent [2] + newax [2];
214     npu = cylPhotonPartition2(emap, mp << 1, newct, newax, d2);
215    
216     /* return total */
217     return npl + npu;
218     }
219    
220    
221    
222     static void cylPhotonPartition (EmissionMap* emap)
223     /* Partition cylindrical source for photon emission */
224     {
225     double d2;
226    
227 greg 2.2 memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
228 greg 2.1 d2 = srcsizerat * thescene.cusize;
229     d2 = PI * emap -> src -> ss2 / (2 * emap -> src -> srad * sqr(d2));
230     d2 *= d2 * DOT(emap -> src -> ss [SU], emap -> src -> ss [SU]);
231    
232     emap -> partitionCnt = 0;
233     emap -> numPartitions = cylPhotonPartition2(emap, 1, emap -> src -> sloc,
234     emap -> src -> ss [SU], d2);
235     emap -> partitionCnt = 0;
236     emap -> partArea = PI * emap -> src -> ss2 / emap -> numPartitions;
237     }
238    
239    
240    
241     static void flatPhotonOrigin (EmissionMap* emap)
242     /* Init emission map with photon origin and associated surface axes on
243     flat light source surface. Also sets source aperture and sampling
244     hemisphere axes at origin */
245     {
246     int i, cent[3], size[3], parr[2];
247     FVECT vpos;
248    
249     cent [0] = cent [1] = cent [2] = 0;
250     size [0] = size [1] = size [2] = emap -> maxPartitions;
251     parr [0] = 0;
252     parr [1] = emap -> partitionCnt;
253    
254     if (!skipparts(cent, size, parr, emap -> partitions))
255     error(CONSISTENCY, "bad source partition in flatPhotonOrigin");
256    
257     vpos [0] = (1 - 2 * pmapRandom(partState)) * size [0] /
258     emap -> maxPartitions;
259     vpos [1] = (1 - 2 * pmapRandom(partState)) * size [1] /
260     emap -> maxPartitions;
261     vpos [2] = 0;
262    
263     for (i = 0; i < 3; i++)
264     vpos [i] += (double)cent [i] / emap -> maxPartitions;
265    
266     /* Get origin */
267     for (i = 0; i < 3; i++)
268     emap -> photonOrg [i] = emap -> src -> sloc [i] +
269     vpos [SU] * emap -> src -> ss [SU][i] +
270     vpos [SV] * emap -> src -> ss [SV][i] +
271     vpos [SW] * emap -> src -> ss [SW][i];
272    
273     /* Get surface axes */
274     VCOPY(emap -> us, emap -> src -> ss [SU]);
275     normalize(emap -> us);
276     VCOPY(emap -> ws, emap -> src -> ss [SW]);
277    
278     if (emap -> port)
279     /* Acts as a photon port; reverse normal as it points INSIDE per
280     * mkillum convention */
281     for (i = 0; i < 3; i++)
282     emap -> ws [i] = -emap -> ws [i];
283    
284     fcross(emap -> vs, emap -> ws, emap -> us);
285    
286     /* Get hemisphere axes & aperture */
287     if (emap -> src -> sflags & SSPOT) {
288     VCOPY(emap -> wh, emap -> src -> sl.s -> aim);
289     i = 0;
290    
291     do {
292     emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
293     emap -> vh [i++] = 1;
294     fcross(emap -> uh, emap -> vh, emap -> wh);
295     } while (normalize(emap -> uh) < FTINY);
296    
297     fcross(emap -> vh, emap -> wh, emap -> uh);
298     emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI);
299     }
300    
301     else {
302     VCOPY(emap -> uh, emap -> us);
303     VCOPY(emap -> vh, emap -> vs);
304     VCOPY(emap -> wh, emap -> ws);
305     emap -> cosThetaMax = 0;
306     }
307     }
308    
309    
310    
311     static void spherePhotonOrigin (EmissionMap* emap)
312     /* Init emission map with photon origin and associated surface axes on
313     spherical light source. Also sets source aperture and sampling
314     hemisphere axes at origin */
315     {
316     int i = 0;
317     unsigned numTheta, numPhi, t, p;
318     RREAL cosTheta, sinTheta, phi;
319    
320     /* Get current partition */
321     numTheta = max(sqrt(2 * emap -> numPartitions / PI) + 0.5, 1);
322     numPhi = 0.5 * PI * numTheta + 0.5;
323    
324     t = emap -> partitionCnt / numPhi;
325     p = emap -> partitionCnt - t * numPhi;
326    
327     emap -> ws [2] = cosTheta = 1 - 2 * (t + pmapRandom(partState)) / numTheta;
328     sinTheta = sqrt(1 - sqr(cosTheta));
329     phi = 2 * PI * (p + pmapRandom(partState)) / numPhi;
330     emap -> ws [0] = cos(phi) * sinTheta;
331     emap -> ws [1] = sin(phi) * sinTheta;
332    
333     if (emap -> port)
334     /* Acts as a photon port; reverse normal as it points INSIDE per
335     * mkillum convention */
336     for (i = 0; i < 3; i++)
337     emap -> ws [i] = -emap -> ws [i];
338    
339     /* Get surface axes us & vs perpendicular to ws */
340     do {
341     emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0;
342     emap -> vs [i++] = 1;
343     fcross(emap -> us, emap -> vs, emap -> ws);
344     } while (normalize(emap -> us) < FTINY);
345    
346     fcross(emap -> vs, emap -> ws, emap -> us);
347    
348     /* Get origin */
349     for (i = 0; i < 3; i++)
350     emap -> photonOrg [i] = emap -> src -> sloc [i] +
351     emap -> src -> srad * emap -> ws [i];
352    
353     /* Get hemisphere axes & aperture */
354     if (emap -> src -> sflags & SSPOT) {
355     VCOPY(emap -> wh, emap -> src -> sl.s -> aim);
356     i = 0;
357    
358     do {
359     emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
360     emap -> vh [i++] = 1;
361     fcross(emap -> uh, emap -> vh, emap -> wh);
362     } while (normalize(emap -> uh) < FTINY);
363    
364     fcross(emap -> vh, emap -> wh, emap -> uh);
365     emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI);
366     }
367    
368     else {
369     VCOPY(emap -> uh, emap -> us);
370     VCOPY(emap -> vh, emap -> vs);
371     VCOPY(emap -> wh, emap -> ws);
372     emap -> cosThetaMax = 0;
373     }
374     }
375    
376    
377    
378     static void sourcePhotonOrigin (EmissionMap* emap)
379     /* Init emission map with photon origin and associated surface axes
380     on scene cube face for distant light source. Also sets source
381     aperture (solid angle) and sampling hemisphere axes at origin */
382     {
383     unsigned long i, partsPerDim, partsPerFace;
384     unsigned face;
385     RREAL du, dv;
386    
387     if (emap -> port) {
388     /* Get origin on photon port */
389     SRCREC *src = emap -> src;
390     emap -> src = emap -> port;
391     photonOrigin [emap -> src -> so -> otype] (emap);
392     emap -> src = src;
393     }
394    
395     else {
396     /* No ports defined, so get origin on scene cube face and SUFFA! */
397     /* Get current face from partition number */
398     partsPerDim = 1 / srcsizerat;
399     partsPerFace = sqr(partsPerDim);
400     face = emap -> partitionCnt / partsPerFace;
401    
402     if (!(emap -> partitionCnt % partsPerFace)) {
403     /* Skipped to a new face; get its normal */
404     emap -> ws [0] = emap -> ws [1] = emap -> ws [2] = 0;
405     emap -> ws [face >> 1] = face & 1 ? 1 : -1;
406    
407     /* Get surface axes us & vs perpendicular to ws */
408     face >>= 1;
409     emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0;
410     emap -> vs [(face + (emap -> ws [face] > 0 ? 2 : 1)) % 3] = 1;
411     fcross(emap -> us, emap -> vs, emap -> ws);
412     }
413    
414     /* Get jittered offsets within face from partition number
415     (in range [-0.5, 0.5]) */
416     i = emap -> partitionCnt % partsPerFace;
417     du = (i / partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5;
418     dv = (i % partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5;
419    
420     /* Jittered destination point within partition */
421     for (i = 0; i < 3; i++)
422     emap -> photonOrg [i] = thescene.cuorg [i] +
423     thescene.cusize * (0.5 + du * emap -> us [i] +
424     dv * emap -> vs [i] +
425     0.5 * emap -> ws [i]);
426     }
427    
428     /* Get hemisphere axes & aperture */
429     VCOPY(emap -> wh, emap -> src -> sloc);
430     i = 0;
431    
432     do {
433     emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
434     emap -> vh [i++] = 1;
435     fcross(emap -> uh, emap -> vh, emap -> wh);
436     } while (normalize(emap -> uh) < FTINY);
437    
438     fcross(emap -> vh, emap -> wh, emap -> uh);
439    
440     /* Get aperture */
441     emap -> cosThetaMax = 1 - emap -> src -> ss2 / (2 * PI);
442     emap -> cosThetaMax = min(1, max(-1, emap -> cosThetaMax));
443     }
444    
445    
446    
447     static void cylPhotonOrigin (EmissionMap* emap)
448     /* Init emission map with photon origin and associated surface axes
449     on cylindrical light source surface. Also sets source aperture
450     and sampling hemisphere axes at origin */
451     {
452     int i, cent[3], size[3], parr[2];
453     FVECT v;
454    
455     cent [0] = cent [1] = cent [2] = 0;
456     size [0] = size [1] = size [2] = emap -> maxPartitions;
457     parr [0] = 0;
458     parr [1] = emap -> partitionCnt;
459    
460     if (!skipparts(cent, size, parr, emap -> partitions))
461     error(CONSISTENCY, "bad source partition in cylPhotonOrigin");
462    
463     v [SU] = 0;
464     v [SV] = (1 - 2 * pmapRandom(partState)) * (double)size [1] /
465     emap -> maxPartitions;
466     v [SW] = (1 - 2 * pmapRandom(partState)) * (double)size [2] /
467     emap -> maxPartitions;
468     normalize(v);
469     v [SU] = (1 - 2 * pmapRandom(partState)) * (double)size [1] /
470     emap -> maxPartitions;
471    
472     for (i = 0; i < 3; i++)
473     v [i] += (double)cent [i] / emap -> maxPartitions;
474    
475     /* Get surface axes */
476     for (i = 0; i < 3; i++)
477     emap -> photonOrg [i] = emap -> ws [i] =
478     (v [SV] * emap -> src -> ss [SV][i] +
479     v [SW] * emap -> src -> ss [SW][i]) / 0.8559;
480    
481     if (emap -> port)
482     /* Acts as a photon port; reverse normal as it points INSIDE per
483     * mkillum convention */
484     for (i = 0; i < 3; i++)
485     emap -> ws [i] = -emap -> ws [i];
486    
487     normalize(emap -> ws);
488     VCOPY(emap -> us, emap -> src -> ss [SU]);
489     normalize(emap -> us);
490     fcross(emap -> vs, emap -> ws, emap -> us);
491    
492     /* Get origin */
493     for (i = 0; i < 3; i++)
494     emap -> photonOrg [i] += v [SU] * emap -> src -> ss [SU][i] +
495     emap -> src -> sloc [i];
496    
497     /* Get hemisphere axes & aperture */
498     if (emap -> src -> sflags & SSPOT) {
499     VCOPY(emap -> wh, emap -> src -> sl.s -> aim);
500     i = 0;
501    
502     do {
503     emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
504     emap -> vh [i++] = 1;
505     fcross(emap -> uh, emap -> vh, emap -> wh);
506     } while (normalize(emap -> uh) < FTINY);
507    
508     fcross(emap -> vh, emap -> wh, emap -> uh);
509     emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI);
510     }
511    
512     else {
513     VCOPY(emap -> uh, emap -> us);
514     VCOPY(emap -> vh, emap -> vs);
515     VCOPY(emap -> wh, emap -> ws);
516     emap -> cosThetaMax = 0;
517     }
518     }
519    
520    
521    
522     void getPhotonPorts ()
523     /* Find geometry declared as photon ports */
524     {
525     OBJECT i;
526     OBJREC* obj;
527    
528     /* Check for missing port modifiers */
529     if (!ambset [0])
530     error(USER, "no photon ports found");
531    
532     for (i = 0; i < nobjects; i++) {
533     obj = objptr(i);
534    
535     if (inset(ambset, obj -> omod)) {
536     /* Add photon port */
537     photonPorts = (SRCREC*)realloc(photonPorts,
538     (numPhotonPorts + 1) *
539     sizeof(SRCREC));
540     if (!photonPorts)
541     error(USER, "can't allocate photon ports");
542    
543     photonPorts [numPhotonPorts].so = obj;
544     photonPorts [numPhotonPorts].sflags = 0;
545    
546     if (!sfun [obj -> otype].of || !sfun[obj -> otype].of -> setsrc)
547     objerror(obj, USER, "illegal photon port");
548    
549     setsource(photonPorts + numPhotonPorts, obj);
550     numPhotonPorts++;
551     }
552     }
553     }
554    
555    
556    
557     static void defaultEmissionFunc (EmissionMap* emap)
558     /* Default behaviour when no emission funcs defined for this source type */
559     {
560     objerror(emap -> src -> so, INTERNAL,
561     "undefined photon emission function");
562     }
563    
564    
565    
566     void initPhotonEmissionFuncs ()
567     /* Init photonPartition[] and photonOrigin[] dispatch tables */
568     {
569     int i;
570    
571     for (i = 0; i < NUMOTYPE; i++)
572     photonPartition [i] = photonOrigin [i] = defaultEmissionFunc;
573    
574     photonPartition [OBJ_FACE] = photonPartition [OBJ_RING] = flatPhotonPartition;
575     photonPartition [OBJ_SOURCE] = sourcePhotonPartition;
576     photonPartition [OBJ_SPHERE] = spherePhotonPartition;
577     photonPartition [OBJ_CYLINDER] = cylPhotonPartition;
578     photonOrigin [OBJ_FACE] = photonOrigin [OBJ_RING] = flatPhotonOrigin;
579     photonOrigin [OBJ_SOURCE] = sourcePhotonOrigin;
580     photonOrigin [OBJ_SPHERE] = spherePhotonOrigin;
581     photonOrigin [OBJ_CYLINDER] = cylPhotonOrigin;
582     }
583    
584    
585    
586     void initPhotonEmission (EmissionMap *emap, float numPdfSamples)
587     /* Initialize photon emission from partitioned light source emap -> src;
588     * this involves integrating the flux emitted from the current photon
589     * origin emap -> photonOrg and setting up a PDF to sample the emission
590     * distribution with numPdfSamples samples */
591     {
592     unsigned i, t, p;
593     double phi, cosTheta, sinTheta, du, dv, dOmega, thetaScale;
594     EmissionSample* sample;
595     const OBJREC* mod = objptr(emap -> src -> so -> omod);
596     static RAY r;
597     #if 0
598     static double lastCosNorm = FHUGE;
599     static SRCREC *lastSrc = NULL, *lastPort = NULL;
600     #endif
601    
602     photonOrigin [emap -> src -> so -> otype] (emap);
603     cosTheta = DOT(emap -> ws, emap -> wh);
604    
605     #if 0
606     if (emap -> src == lastSrc && emap -> port == lastPort &&
607     (emap -> src -> sflags & SDISTANT || mod -> omod == OVOID) &&
608     cosTheta == lastCosNorm)
609     /* Same source, port, and aperture-normal angle, and source is
610     either distant (and thus translationally invariant) or has
611     no modifier --> flux unchanged */
612     /* BUG: this optimisation ignores partial occlusion of ports and
613     can lead to erroneous "zero emission" bailouts.
614     It can also lead to bias with modifiers exhibiting high variance!
615     Disabled for now -- RS 12/13 */
616     return;
617    
618     lastSrc = emap -> src;
619     lastPort = emap -> port;
620     lastCosNorm = cosTheta;
621     #endif
622    
623     /* Need to recompute flux & PDF */
624     setcolor(emap -> partFlux, 0, 0, 0);
625     emap -> cdf = 0;
626     emap -> numSamples = 0;
627    
628     if (cosTheta <= 0 &&
629     sqrt(1 - sqr(cosTheta)) <= emap -> cosThetaMax + FTINY)
630     /* Aperture below surface; no emission from current origin */
631     return;
632    
633     if (mod -> omod == OVOID && !emap -> port &&
634     (cosTheta >= 1 - FTINY || (emap -> src -> sflags & SDISTANT &&
635     acos(cosTheta) + acos(emap -> cosThetaMax) <= 0.5 * PI))) {
636     /* Source is unmodified and has no port (which requires testing for
637     occlusion), and is either local with normal aligned aperture or
638     distant with aperture above surface; analytical flux & PDF */
639     setcolor(emap -> partFlux, mod -> oargs.farg [0],
640     mod -> oargs.farg [1], mod -> oargs.farg [2]);
641    
642     /* Multiply radiance by Omega * dA to get flux */
643     scalecolor(emap -> partFlux,
644     PI * cosTheta * (1 - sqr(max(emap -> cosThetaMax, 0))) *
645     emap -> partArea);
646     }
647    
648     else {
649     /* Source is either modified, has a port, is local with off-normal
650     aperture, or distant with aperture partly below surface; get flux
651     via numerical integration */
652     thetaScale = (1 - emap -> cosThetaMax);
653    
654     /* Figga out numba of aperture samples for integration;
655     numTheta / numPhi ratio is 1 / PI */
656     du = sqrt(pdfSamples * 2 * thetaScale);
657     emap -> numTheta = max(du + 0.5, 1);
658     emap -> numPhi = max(PI * du + 0.5, 1);
659    
660     dOmega = 2 * PI * thetaScale / (emap -> numTheta * emap -> numPhi);
661     thetaScale /= emap -> numTheta;
662    
663     /* Allocate PDF, baby */
664     sample = emap -> samples = (EmissionSample*)
665     realloc(emap -> samples,
666     sizeof(EmissionSample) *
667     emap -> numTheta * emap -> numPhi);
668     if (!emap -> samples)
669     error(USER, "can't allocate emission PDF");
670    
671     VCOPY(r.rorg, emap -> photonOrg);
672     VCOPY(r.rop, emap -> photonOrg);
673     r.rmax = FHUGE;
674    
675     for (t = 0; t < emap -> numTheta; t++) {
676     for (p = 0; p < emap -> numPhi; p++) {
677     /* This uniform mapping handles 0 <= thetaMax <= PI */
678     cosTheta = 1 - (t + pmapRandom(emitState)) * thetaScale;
679     sinTheta = sqrt(1 - sqr(cosTheta));
680     phi = 2 * PI * (p + pmapRandom(emitState)) / emap -> numPhi;
681     du = cos(phi) * sinTheta;
682     dv = sin(phi) * sinTheta;
683     rayorigin(&r, PRIMARY, NULL, NULL);
684    
685     for (i = 0; i < 3; i++)
686     r.rdir [i] = du * emap -> uh [i] + dv * emap -> vh [i] +
687     cosTheta * emap -> wh [i];
688    
689     /* Sample behind surface? */
690     VCOPY(r.ron, emap -> ws);
691     if ((r.rod = DOT(r.rdir, r.ron)) <= 0)
692     continue;
693    
694     /* Get radiance emitted in this direction; to get flux we
695     multiply by cos(theta_surface), dOmega, and dA. Ray
696     is directed towards light source for raytexture(). */
697     if (!(emap -> src -> sflags & SDISTANT))
698     for (i = 0; i < 3; i++)
699     r.rdir [i] = -r.rdir [i];
700    
701     /* Port occluded in this direction? */
702     if (emap -> port && localhit(&r, &thescene))
703     continue;
704    
705     raytexture(&r, mod -> omod);
706     setcolor(r.rcol, mod -> oargs.farg [0], mod -> oargs.farg [1],
707     mod -> oargs.farg [2]);
708     multcolor(r.rcol, r.pcol);
709    
710     /* Multiply by cos(theta_surface) */
711     scalecolor(r.rcol, r.rod);
712    
713     /* Add PDF sample if nonzero; importance info for photon emission
714     * could go here... */
715     if (colorAvg(r.rcol)) {
716     copycolor(sample -> pdf, r.rcol);
717     sample -> cdf = emap -> cdf += colorAvg(sample -> pdf);
718     sample -> theta = t;
719     sample++ -> phi = p;
720     emap -> numSamples++;
721     addcolor(emap -> partFlux, r.rcol);
722     }
723     }
724     }
725    
726     /* Multiply by dOmega * dA */
727     scalecolor(emap -> partFlux, dOmega * emap -> partArea);
728     }
729     }
730    
731    
732    
733     #define vomitPhoton emitPhoton
734     #define bluarrrghPhoton vomitPhoton
735    
736     void emitPhoton (const EmissionMap* emap, RAY* ray)
737     /* Emit photon from current partition emap -> partitionCnt based on
738     emission distribution. Returns new photon ray. */
739     {
740     unsigned long i, lo, hi;
741     const EmissionSample* sample = emap -> samples;
742     RREAL du, dv, cosTheta, cosThetaSqr, sinTheta, phi;
743     const OBJREC* mod = objptr(emap -> src -> so -> omod);
744    
745     /* Choose a new origin within current partition for every
746     emitted photon to break up clustering artifacts */
747     photonOrigin [emap -> src -> so -> otype] ((EmissionMap*)emap);
748     rayorigin(ray, PRIMARY, NULL, NULL);
749     ray -> rmax = FHUGE;
750    
751     if (!emap -> numSamples) {
752     /* Source is unmodified and has no port, and either local with
753     normal aligned aperture, or distant with aperture above surface;
754     use cosine weighted distribution */
755     cosThetaSqr = 1 - pmapRandom(emitState) *
756     (1 - sqr(max(emap -> cosThetaMax, 0)));
757     cosTheta = sqrt(cosThetaSqr);
758     sinTheta = sqrt(1 - cosThetaSqr);
759     phi = 2 * PI * pmapRandom(emitState);
760     setcolor(ray -> rcol, mod -> oargs.farg [0], mod -> oargs.farg [1],
761     mod -> oargs.farg [2]);
762     }
763    
764     else {
765     /* Source is either modified, has a port, is local with off-normal
766     aperture, or distant with aperture partly below surface; choose
767     direction from constructed cumulative distribution function with
768     Monte Carlo inversion using binary search. */
769     du = pmapRandom(emitState) * emap -> cdf;
770     lo = 1;
771     hi = emap -> numSamples;
772    
773     while (hi > lo) {
774     i = (lo + hi) >> 1;
775     sample = emap -> samples + i - 1;
776    
777     if (sample -> cdf >= du)
778     hi = i;
779     if (sample -> cdf < du)
780     lo = i + 1;
781     }
782    
783     /* This is a uniform mapping, mon */
784     cosTheta = 1 - (sample -> theta + pmapRandom(emitState)) *
785     (1 - emap -> cosThetaMax) / emap -> numTheta;
786     sinTheta = sqrt(1 - sqr(cosTheta));
787     phi = 2 * PI * (sample -> phi + pmapRandom(emitState)) / emap -> numPhi;
788     copycolor(ray -> rcol, sample -> pdf);
789     }
790    
791     /* Normalize photon flux so that average over RGB is 1 */
792     colorNorm(ray -> rcol);
793    
794     VCOPY(ray -> rorg, emap -> photonOrg);
795     du = cos(phi) * sinTheta;
796     dv = sin(phi) * sinTheta;
797    
798     for (i = 0; i < 3; i++)
799     ray -> rdir [i] = du * emap -> uh [i] + dv * emap -> vh [i] +
800     cosTheta * emap -> wh [i];
801    
802     if (emap -> src -> sflags & SDISTANT)
803     /* Distant source; reverse ray direction to point into the scene. */
804     for (i = 0; i < 3; i++)
805     ray -> rdir [i] = -ray -> rdir [i];
806    
807     if (emap -> port)
808     /* Photon emitted from port; move origin just behind port so it
809     will be scattered */
810     for (i = 0; i < 3; i++)
811     ray -> rorg [i] -= 2 * FTINY * ray -> rdir [i];
812    
813     /* Assign emitting light source index */
814     ray -> rsrc = emap -> src - source;
815     }
816    
817    
818    
819     void multDirectPmap (RAY *r)
820     /* Factor irradiance from direct photons into r -> rcol; interface to
821     * direct() */
822     {
823     COLOR photonIrrad;
824    
825     /* Lookup direct photon irradiance */
826     (directPmap -> lookup)(directPmap, r, photonIrrad);
827    
828     /* Multiply with coefficient in ray */
829     multcolor(r -> rcol, photonIrrad);
830    
831     return;
832     }
833    
834    
835    
836     void inscatterVolumePmap (RAY *r, COLOR inscatter)
837     /* Add inscattering from volume photon map; interface to srcscatter() */
838     {
839     /* Add ambient in-scattering via lookup callback */
840     (volumePmap -> lookup)(volumePmap, r, inscatter);
841     }