ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapsrc.c
Revision: 2.7
Committed: Tue Aug 18 18:45:55 2015 UTC (8 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +4 -1 lines
Log Message:
Added missing RCSid forgotten during initial check-in

File Contents

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