ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapsrc.c
Revision: 2.1
Committed: Tue Feb 24 19:39:27 2015 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial check-in of photon map addition by Roland Schregle

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