ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapsrc.c
Revision: 2.10
Committed: Tue Sep 29 18:16:34 2015 UTC (8 years, 7 months ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.9: +5 -4 lines
Log Message:
Improved handling of mixtures and photon ports with photon mapping

File Contents

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