ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapsrc.c
Revision: 2.2
Committed: Tue Apr 21 19:16:51 2015 UTC (9 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +9 -8 lines
Log Message:
Fixes for Windows and photon map

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