ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapsrc.c
Revision: 2.6
Committed: Thu May 28 12:27:22 2015 UTC (8 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +2 -2 lines
Log Message:
Bug fix for uninitialized sources set to -1

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 greg 2.6 $Id: pmapsrc.c,v 2.5 2015/05/22 11:38:26 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 greg 2.4 const OBJREC* mod = findmaterial(emap -> src -> so);
596 greg 2.1 static RAY r;
597     #if 0
598     static double lastCosNorm = FHUGE;
599     static SRCREC *lastSrc = NULL, *lastPort = NULL;
600     #endif
601    
602 greg 2.4 setcolor(emap -> partFlux, 0, 0, 0);
603    
604 greg 2.1 photonOrigin [emap -> src -> so -> otype] (emap);
605     cosTheta = DOT(emap -> ws, emap -> wh);
606    
607     #if 0
608     if (emap -> src == lastSrc && emap -> port == lastPort &&
609     (emap -> src -> sflags & SDISTANT || mod -> omod == OVOID) &&
610     cosTheta == lastCosNorm)
611     /* Same source, port, and aperture-normal angle, and source is
612     either distant (and thus translationally invariant) or has
613     no modifier --> flux unchanged */
614     /* BUG: this optimisation ignores partial occlusion of ports and
615     can lead to erroneous "zero emission" bailouts.
616     It can also lead to bias with modifiers exhibiting high variance!
617     Disabled for now -- RS 12/13 */
618     return;
619    
620     lastSrc = emap -> src;
621     lastPort = emap -> port;
622     lastCosNorm = cosTheta;
623     #endif
624    
625     /* Need to recompute flux & PDF */
626     emap -> cdf = 0;
627     emap -> numSamples = 0;
628    
629     if (cosTheta <= 0 &&
630     sqrt(1 - sqr(cosTheta)) <= emap -> cosThetaMax + FTINY)
631     /* Aperture below surface; no emission from current origin */
632     return;
633    
634     if (mod -> omod == OVOID && !emap -> port &&
635     (cosTheta >= 1 - FTINY || (emap -> src -> sflags & SDISTANT &&
636     acos(cosTheta) + acos(emap -> cosThetaMax) <= 0.5 * PI))) {
637     /* Source is unmodified and has no port (which requires testing for
638     occlusion), and is either local with normal aligned aperture or
639     distant with aperture above surface; analytical flux & PDF */
640     setcolor(emap -> partFlux, mod -> oargs.farg [0],
641     mod -> oargs.farg [1], mod -> oargs.farg [2]);
642    
643     /* Multiply radiance by Omega * dA to get flux */
644     scalecolor(emap -> partFlux,
645     PI * cosTheta * (1 - sqr(max(emap -> cosThetaMax, 0))) *
646     emap -> partArea);
647     }
648    
649     else {
650     /* Source is either modified, has a port, is local with off-normal
651     aperture, or distant with aperture partly below surface; get flux
652     via numerical integration */
653     thetaScale = (1 - emap -> cosThetaMax);
654    
655     /* Figga out numba of aperture samples for integration;
656     numTheta / numPhi ratio is 1 / PI */
657     du = sqrt(pdfSamples * 2 * thetaScale);
658     emap -> numTheta = max(du + 0.5, 1);
659     emap -> numPhi = max(PI * du + 0.5, 1);
660    
661     dOmega = 2 * PI * thetaScale / (emap -> numTheta * emap -> numPhi);
662     thetaScale /= emap -> numTheta;
663    
664     /* Allocate PDF, baby */
665     sample = emap -> samples = (EmissionSample*)
666     realloc(emap -> samples,
667     sizeof(EmissionSample) *
668     emap -> numTheta * emap -> numPhi);
669     if (!emap -> samples)
670     error(USER, "can't allocate emission PDF");
671    
672     VCOPY(r.rorg, emap -> photonOrg);
673     VCOPY(r.rop, emap -> photonOrg);
674 greg 2.5 r.rmax = 0;
675 greg 2.1
676     for (t = 0; t < emap -> numTheta; t++) {
677     for (p = 0; p < emap -> numPhi; p++) {
678     /* This uniform mapping handles 0 <= thetaMax <= PI */
679     cosTheta = 1 - (t + pmapRandom(emitState)) * thetaScale;
680     sinTheta = sqrt(1 - sqr(cosTheta));
681     phi = 2 * PI * (p + pmapRandom(emitState)) / emap -> numPhi;
682     du = cos(phi) * sinTheta;
683     dv = sin(phi) * sinTheta;
684     rayorigin(&r, PRIMARY, NULL, NULL);
685    
686     for (i = 0; i < 3; i++)
687     r.rdir [i] = du * emap -> uh [i] + dv * emap -> vh [i] +
688     cosTheta * emap -> wh [i];
689    
690     /* Sample behind surface? */
691     VCOPY(r.ron, emap -> ws);
692     if ((r.rod = DOT(r.rdir, r.ron)) <= 0)
693     continue;
694    
695     /* Get radiance emitted in this direction; to get flux we
696     multiply by cos(theta_surface), dOmega, and dA. Ray
697     is directed towards light source for raytexture(). */
698     if (!(emap -> src -> sflags & SDISTANT))
699     for (i = 0; i < 3; i++)
700     r.rdir [i] = -r.rdir [i];
701    
702     /* Port occluded in this direction? */
703     if (emap -> port && localhit(&r, &thescene))
704     continue;
705    
706     raytexture(&r, mod -> omod);
707     setcolor(r.rcol, mod -> oargs.farg [0], mod -> oargs.farg [1],
708     mod -> oargs.farg [2]);
709     multcolor(r.rcol, r.pcol);
710    
711     /* Multiply by cos(theta_surface) */
712     scalecolor(r.rcol, r.rod);
713    
714     /* Add PDF sample if nonzero; importance info for photon emission
715     * could go here... */
716     if (colorAvg(r.rcol)) {
717     copycolor(sample -> pdf, r.rcol);
718     sample -> cdf = emap -> cdf += colorAvg(sample -> pdf);
719     sample -> theta = t;
720     sample++ -> phi = p;
721     emap -> numSamples++;
722     addcolor(emap -> partFlux, r.rcol);
723     }
724     }
725     }
726    
727     /* Multiply by dOmega * dA */
728     scalecolor(emap -> partFlux, dOmega * emap -> partArea);
729     }
730     }
731    
732    
733    
734     #define vomitPhoton emitPhoton
735     #define bluarrrghPhoton vomitPhoton
736    
737     void emitPhoton (const EmissionMap* emap, RAY* ray)
738     /* Emit photon from current partition emap -> partitionCnt based on
739     emission distribution. Returns new photon ray. */
740     {
741     unsigned long i, lo, hi;
742     const EmissionSample* sample = emap -> samples;
743     RREAL du, dv, cosTheta, cosThetaSqr, sinTheta, phi;
744 greg 2.4 const OBJREC* mod = findmaterial(emap -> src -> so);
745 greg 2.1
746     /* Choose a new origin within current partition for every
747     emitted photon to break up clustering artifacts */
748     photonOrigin [emap -> src -> so -> otype] ((EmissionMap*)emap);
749 greg 2.5 /* If we have a local glow source with a maximum radius, then
750     restrict our photon to the specified distance (otherwise no limit) */
751 greg 2.6 if (mod -> otype == MAT_GLOW && !(emap -> src -> sflags & SDISTANT)
752 greg 2.5 && mod -> oargs.farg[3] > FTINY)
753     ray -> rmax = mod -> oargs.farg[3];
754     else
755     ray -> rmax = 0;
756 greg 2.1 rayorigin(ray, PRIMARY, NULL, NULL);
757    
758     if (!emap -> numSamples) {
759     /* Source is unmodified and has no port, and either local with
760     normal aligned aperture, or distant with aperture above surface;
761     use cosine weighted distribution */
762     cosThetaSqr = 1 - pmapRandom(emitState) *
763     (1 - sqr(max(emap -> cosThetaMax, 0)));
764     cosTheta = sqrt(cosThetaSqr);
765     sinTheta = sqrt(1 - cosThetaSqr);
766     phi = 2 * PI * pmapRandom(emitState);
767     setcolor(ray -> rcol, mod -> oargs.farg [0], mod -> oargs.farg [1],
768     mod -> oargs.farg [2]);
769     }
770    
771     else {
772     /* Source is either modified, has a port, is local with off-normal
773     aperture, or distant with aperture partly below surface; choose
774     direction from constructed cumulative distribution function with
775     Monte Carlo inversion using binary search. */
776     du = pmapRandom(emitState) * emap -> cdf;
777     lo = 1;
778     hi = emap -> numSamples;
779    
780     while (hi > lo) {
781     i = (lo + hi) >> 1;
782     sample = emap -> samples + i - 1;
783    
784     if (sample -> cdf >= du)
785     hi = i;
786     if (sample -> cdf < du)
787     lo = i + 1;
788     }
789    
790     /* This is a uniform mapping, mon */
791     cosTheta = 1 - (sample -> theta + pmapRandom(emitState)) *
792     (1 - emap -> cosThetaMax) / emap -> numTheta;
793     sinTheta = sqrt(1 - sqr(cosTheta));
794     phi = 2 * PI * (sample -> phi + pmapRandom(emitState)) / emap -> numPhi;
795     copycolor(ray -> rcol, sample -> pdf);
796     }
797    
798     /* Normalize photon flux so that average over RGB is 1 */
799     colorNorm(ray -> rcol);
800    
801     VCOPY(ray -> rorg, emap -> photonOrg);
802     du = cos(phi) * sinTheta;
803     dv = sin(phi) * sinTheta;
804    
805     for (i = 0; i < 3; i++)
806     ray -> rdir [i] = du * emap -> uh [i] + dv * emap -> vh [i] +
807     cosTheta * emap -> wh [i];
808    
809     if (emap -> src -> sflags & SDISTANT)
810     /* Distant source; reverse ray direction to point into the scene. */
811     for (i = 0; i < 3; i++)
812     ray -> rdir [i] = -ray -> rdir [i];
813    
814     if (emap -> port)
815     /* Photon emitted from port; move origin just behind port so it
816     will be scattered */
817     for (i = 0; i < 3; i++)
818     ray -> rorg [i] -= 2 * FTINY * ray -> rdir [i];
819    
820     /* Assign emitting light source index */
821     ray -> rsrc = emap -> src - source;
822     }
823    
824    
825    
826     void multDirectPmap (RAY *r)
827     /* Factor irradiance from direct photons into r -> rcol; interface to
828     * direct() */
829     {
830     COLOR photonIrrad;
831    
832     /* Lookup direct photon irradiance */
833     (directPmap -> lookup)(directPmap, r, photonIrrad);
834    
835     /* Multiply with coefficient in ray */
836     multcolor(r -> rcol, photonIrrad);
837    
838     return;
839     }
840    
841    
842    
843     void inscatterVolumePmap (RAY *r, COLOR inscatter)
844     /* Add inscattering from volume photon map; interface to srcscatter() */
845     {
846     /* Add ambient in-scattering via lookup callback */
847     (volumePmap -> lookup)(volumePmap, r, inscatter);
848     }