ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapsrc.c
Revision: 2.9
Committed: Tue Sep 22 15:08:31 2015 UTC (9 years, 7 months ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.8: +5 -2 lines
Log Message:
Improved photon port handling -- only actual surfaces are accepted

File Contents

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