ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmap.c
Revision: 2.3
Committed: Wed Apr 22 20:28:16 2015 UTC (9 years, 1 month ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.2: +18 -10 lines
Log Message:
Disabled use of SIGCONT if undefined (e.g. on WIN32), removed extra space in
rpict reports after biascomp stats

File Contents

# User Rev Content
1 greg 2.1 /*
2     ==================================================================
3     Photon map main module
4    
5     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
6     (c) Fraunhofer Institute for Solar Energy Systems,
7     Lucerne University of Applied Sciences & Arts
8     ==================================================================
9    
10 rschregle 2.3 $Id: pmap.c,v 2.2 2015/04/21 19:16:51 greg Exp $
11 greg 2.1 */
12    
13    
14    
15     #include "pmap.h"
16     #include "pmapmat.h"
17     #include "pmapsrc.h"
18     #include "pmaprand.h"
19     #include "pmapio.h"
20     #include "pmapbias.h"
21     #include "pmapdiag.h"
22     #include "otypes.h"
23     #include <time.h>
24     #include <sys/stat.h>
25    
26    
27    
28     extern char *octname;
29    
30 rschregle 2.3 static char PmapRevision [] = "$Revision: 2.2 $";
31 greg 2.1
32    
33    
34     /* Photon map lookup functions per type */
35     void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
36     photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
37     photonDensity, NULL
38     };
39    
40    
41    
42     void colorNorm (COLOR c)
43     /* Normalise colour channels to average of 1 */
44     {
45     const float avg = colorAvg(c);
46    
47     if (!avg)
48     return;
49    
50     c [0] /= avg;
51     c [1] /= avg;
52     c [2] /= avg;
53     }
54    
55    
56    
57     void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm)
58     {
59     unsigned t;
60     struct stat octstat, pmstat;
61     PhotonMap *pm;
62     PhotonMapType type;
63    
64     for (t = 0; t < NUM_PMAP_TYPES; t++)
65     if (setPmapParam(&pm, parm + t)) {
66     /* Check if photon map newer than octree */
67     if (!stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) &&
68     octstat.st_mtime > pmstat.st_mtime) {
69     sprintf(errmsg, "photon map in file %s may be stale",
70     pm -> fileName);
71     error(USER, errmsg);
72     }
73    
74     /* Load photon map from file and get its type */
75     if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE)
76     error(USER, "failed loading photon map");
77    
78     /* Assign to appropriate photon map type (deleting previously
79     * loaded photon map of same type if necessary) */
80     if (pmaps [type]) {
81     deletePhotons(pmaps [type]);
82     free(pmaps [type]);
83     }
84     pmaps [type] = pm;
85    
86     /* Check for invalid density estimate bandwidth */
87     if (pm -> maxGather > pm -> heapSize) {
88     error(WARNING, "adjusting density estimate bandwidth");
89     pm -> minGather = pm -> maxGather = pm -> heapSize;
90     }
91     }
92     }
93    
94    
95    
96     void savePmaps (const PhotonMap **pmaps, int argc, char **argv)
97     {
98     unsigned t;
99    
100     for (t = 0; t < NUM_PMAP_TYPES; t++) {
101     if (pmaps [t])
102     savePhotonMap(pmaps [t], pmaps [t] -> fileName, t, argc, argv);
103     }
104     }
105    
106    
107    
108     void cleanUpPmaps (PhotonMap **pmaps)
109     {
110     unsigned t;
111    
112     for (t = 0; t < NUM_PMAP_TYPES; t++) {
113     if (pmaps [t]) {
114     deletePhotons(pmaps [t]);
115     free(pmaps [t]);
116     }
117     }
118     }
119    
120    
121    
122     static int photonParticipate (RAY *ray)
123     /* Trace photon through participating medium. Returns 1 if passed through,
124     or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */
125     {
126     int i;
127     RREAL cosTheta, cosPhi, du, dv;
128     const float cext = colorAvg(ray -> cext),
129     albedo = colorAvg(ray -> albedo);
130     FVECT u, v;
131     COLOR cvext;
132    
133     /* Mean free distance until interaction with medium */
134     ray -> rmax = -log(pmapRandom(mediumState)) / cext;
135    
136     while (!localhit(ray, &thescene)) {
137     setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]),
138     exp(-ray -> rmax * ray -> cext [1]),
139     exp(-ray -> rmax * ray -> cext [2]));
140    
141     /* Modify ray color and normalise */
142     multcolor(ray -> rcol, cvext);
143     colorNorm(ray -> rcol);
144     VCOPY(ray -> rorg, ray -> rop);
145    
146     if (albedo > FTINY)
147     /* Add to volume photon map */
148     if (ray -> rlvl > 0) addPhoton(volumePmap, ray);
149    
150     /* Absorbed? */
151     if (pmapRandom(rouletteState) > albedo) return 0;
152    
153     /* Colour bleeding without attenuation (?) */
154     multcolor(ray -> rcol, ray -> albedo);
155     scalecolor(ray -> rcol, 1 / albedo);
156    
157     /* Scatter photon */
158     cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1
159     : 1 / (2 * ray -> gecc) *
160     (1 + ray -> gecc * ray -> gecc -
161     (1 - ray -> gecc * ray -> gecc) /
162     (1 - ray -> gecc + 2 * ray -> gecc *
163     pmapRandom(scatterState)));
164    
165     cosPhi = cos(2 * PI * pmapRandom(scatterState));
166     du = dv = sqrt(1 - cosTheta * cosTheta); /* sin(theta) */
167     du *= cosPhi;
168     dv *= sqrt(1 - cosPhi * cosPhi); /* sin(phi) */
169    
170     /* Get axes u & v perpendicular to photon direction */
171     i = 0;
172     do {
173     v [0] = v [1] = v [2] = 0;
174     v [i++] = 1;
175     fcross(u, v, ray -> rdir);
176     } while (normalize(u) < FTINY);
177     fcross(v, ray -> rdir, u);
178    
179     for (i = 0; i < 3; i++)
180     ray -> rdir [i] = du * u [i] + dv * v [i] +
181     cosTheta * ray -> rdir [i];
182     ray -> rlvl++;
183     ray -> rmax = -log(pmapRandom(mediumState)) / cext;
184     }
185    
186     setcolor(cvext, exp(-ray -> rot * ray -> cext [0]),
187     exp(-ray -> rot * ray -> cext [1]),
188     exp(-ray -> rot * ray -> cext [2]));
189    
190     /* Modify ray color and normalise */
191     multcolor(ray -> rcol, cvext);
192     colorNorm(ray -> rcol);
193    
194     /* Passed through medium */
195     return 1;
196     }
197    
198    
199    
200     void tracePhoton (RAY *ray)
201     /* Follow photon as it bounces around... */
202     {
203     long mod;
204     OBJREC* mat;
205    
206     if (ray -> rlvl > photonMaxBounce) {
207     error(WARNING, "runaway photon!");
208     return;
209     }
210    
211     if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
212     return;
213    
214     if (localhit(ray, &thescene)) {
215     mod = ray -> ro -> omod;
216    
217     if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) {
218     /* Transfer ray if modifier is VOID or clipped within antimatta */
219     RAY tray;
220     photonRay(ray, &tray, PMAP_XFER, NULL);
221     tracePhoton(&tray);
222     }
223     else {
224     /* Scatter for modifier material */
225     mat = objptr(mod);
226     photonScatter [mat -> otype] (mat, ray);
227     }
228     }
229     }
230    
231    
232    
233     static void preComputeGlobal (PhotonMap *pmap)
234     /* Precompute irradiance from global photons for final gathering using
235     the first finalGather * pmap -> heapSize photons in the heap. Returns
236     new heap with precomputed photons. */
237     {
238     unsigned long i, nuHeapSize;
239     unsigned j;
240     Photon *nuHeap, *p;
241     COLOR irrad;
242     RAY ray;
243     float nuMinPos [3], nuMaxPos [3];
244    
245     repComplete = nuHeapSize = finalGather * pmap -> heapSize;
246    
247     if (photonRepTime) {
248     sprintf(errmsg,
249     "Precomputing irradiance for %ld global photons...\n",
250     nuHeapSize);
251     eputs(errmsg);
252     fflush(stderr);
253     }
254    
255     p = nuHeap = (Photon*)malloc(nuHeapSize * sizeof(Photon));
256     if (!nuHeap)
257     error(USER, "can't allocate photon heap");
258    
259     for (j = 0; j <= 2; j++) {
260     nuMinPos [j] = FHUGE;
261     nuMaxPos [j] = -FHUGE;
262     }
263    
264     /* Record start time, baby */
265     repStartTime = time(NULL);
266 rschregle 2.3 #ifdef SIGCONT
267     signal(SIGCONT, pmapPreCompReport);
268     #endif
269 greg 2.1 repProgress = 0;
270 greg 2.2 memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon));
271 greg 2.1
272     for (i = 0, p = nuHeap; i < nuHeapSize; i++, p++) {
273     ray.ro = NULL;
274     VCOPY(ray.rop, p -> pos);
275    
276     /* Update min and max positions & set ray normal */
277     for (j = 0; j < 3; j++) {
278     if (p -> pos [j] < nuMinPos [j]) nuMinPos [j] = p -> pos [j];
279     if (p -> pos [j] > nuMaxPos [j]) nuMaxPos [j] = p -> pos [j];
280     ray.ron [j] = p -> norm [j] / 127.0;
281     }
282    
283     photonDensity(pmap, &ray, irrad);
284     setPhotonFlux(p, irrad);
285     repProgress++;
286    
287     if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
288     pmapPreCompReport();
289 rschregle 2.3 #ifdef SIGCONT
290     else signal(SIGCONT, pmapPreCompReport);
291     #endif
292 greg 2.1 }
293    
294 rschregle 2.3 #ifdef SIGCONT
295     signal(SIGCONT, SIG_DFL);
296     #endif
297 greg 2.1
298     /* Replace & rebuild heap */
299     free(pmap -> heap);
300     pmap -> heap = nuHeap;
301     pmap -> heapSize = pmap -> heapEnd = nuHeapSize;
302     VCOPY(pmap -> minPos, nuMinPos);
303     VCOPY(pmap -> maxPos, nuMaxPos);
304    
305     if (photonRepTime) {
306     eputs("Rebuilding global photon heap...\n");
307     fflush(stderr);
308     }
309    
310     balancePhotons(pmap, NULL);
311     }
312    
313    
314    
315     void distribPhotons (PhotonMap **pmaps)
316     {
317     EmissionMap emap;
318     char errmsg2 [128];
319     unsigned t, srcIdx, passCnt = 0, prePassCnt = 0;
320     double totalFlux = 0;
321     PhotonMap *pm;
322    
323     for (t = 0; t < NUM_PMAP_TYPES && !photonMaps [t]; t++);
324     if (t >= NUM_PMAP_TYPES)
325     error(USER, "no photon maps defined");
326    
327     if (!nsources)
328     error(USER, "no light sources");
329    
330     /* ===================================================================
331     * INITIALISATION - Set up emission and scattering funcs
332     * =================================================================== */
333     emap.samples = NULL;
334     emap.maxPartitions = MAXSPART;
335     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
336     if (!emap.partitions)
337     error(INTERNAL, "can't allocate source partitions");
338    
339     /* Initialise all defined photon maps */
340     for (t = 0; t < NUM_PMAP_TYPES; t++)
341     initPhotonMap(photonMaps [t], t);
342    
343     initPhotonEmissionFuncs();
344     initPhotonScatterFuncs();
345    
346     /* Get photon ports if specified */
347     if (ambincl == 1)
348     getPhotonPorts();
349    
350     /* Get photon sensor modifiers */
351     getPhotonSensors(photonSensorList);
352    
353     /* Seed RNGs for photon distribution */
354     pmapSeed(randSeed, partState);
355     pmapSeed(randSeed, emitState);
356     pmapSeed(randSeed, cntState);
357     pmapSeed(randSeed, mediumState);
358     pmapSeed(randSeed, scatterState);
359     pmapSeed(randSeed, rouletteState);
360    
361     if (photonRepTime)
362     eputs("\n");
363    
364     /* ===================================================================
365     * FLUX INTEGRATION - Get total photon flux from light sources
366     * =================================================================== */
367     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
368     unsigned portCnt = 0;
369     emap.src = source + srcIdx;
370    
371     do {
372     emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
373     : NULL;
374     photonPartition [emap.src -> so -> otype] (&emap);
375    
376     if (photonRepTime) {
377     sprintf(errmsg, "Integrating flux from source %s ",
378     source [srcIdx].so -> oname);
379    
380     if (emap.port) {
381     sprintf(errmsg2, "via port %s ",
382     photonPorts [portCnt].so -> oname);
383     strcat(errmsg, errmsg2);
384     }
385    
386     sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
387     strcat(errmsg, errmsg2);
388     eputs(errmsg);
389     fflush(stderr);
390     }
391    
392     for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
393     emap.partitionCnt++) {
394     initPhotonEmission(&emap, pdfSamples);
395     totalFlux += colorAvg(emap.partFlux);
396     }
397    
398     portCnt++;
399     } while (portCnt < numPhotonPorts);
400     }
401    
402     if (totalFlux < FTINY)
403     error(USER, "zero flux from light sources");
404    
405     /* Record start time and enable progress report signal handler */
406     repStartTime = time(NULL);
407 rschregle 2.3 #ifdef SIGCONT
408     signal(SIGCONT, pmapDistribReport);
409     #endif
410 greg 2.1 repProgress = prePassCnt = 0;
411    
412     if (photonRepTime)
413     eputs("\n");
414    
415     /* ===================================================================
416     * 2-PASS PHOTON DISTRIBUTION
417     * Pass 1 (pre): emit fraction of target photon count
418     * Pass 2 (main): based on outcome of pass 1, estimate remaining number
419     * of photons to emit to approximate target count
420     * =================================================================== */
421     do {
422     double numEmit;
423    
424     if (!passCnt) {
425     /* INIT PASS 1 */
426     /* Skip if no photons contributed after sufficient iterations; make
427     * it clear to user which photon maps are missing so (s)he can
428     * check the scene geometry and materials */
429     if (++prePassCnt > maxPreDistrib) {
430     sprintf(errmsg, "too many prepasses");
431    
432     for (t = 0; t < NUM_PMAP_TYPES; t++)
433     if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
434     sprintf(errmsg2, ", no %s photons stored", pmapName [t]);
435     strcat(errmsg, errmsg2);
436     }
437    
438     error(USER, errmsg);
439     break;
440     }
441    
442     /* Num to emit is fraction of minimum target count */
443     numEmit = FHUGE;
444    
445     for (t = 0; t < NUM_PMAP_TYPES; t++)
446     if (photonMaps [t])
447     numEmit = min(photonMaps [t] -> distribTarget, numEmit);
448    
449     numEmit *= preDistrib;
450     }
451    
452     else {
453     /* INIT PASS 2 */
454     /* Based on the outcome of the predistribution we can now estimate
455     * how many more photons we have to emit for each photon map to
456     * meet its respective target count. This value is clamped to 0 in
457     * case the target has already been exceeded in the pass 1. Note
458     * repProgress is the number of photons emitted thus far, while
459     * heapEnd is the number of photons stored in each photon map. */
460     double maxDistribRatio = 0;
461    
462     /* Set the distribution ratio for each map; this indicates how many
463     * photons of each respective type are stored per emitted photon,
464     * and is used as probability for storing a photon by addPhoton().
465     * Since this biases the photon density, addPhoton() promotes the
466     * flux of stored photons to compensate. */
467     for (t = 0; t < NUM_PMAP_TYPES; t++)
468     if ((pm = photonMaps [t])) {
469     pm -> distribRatio = (double)pm -> distribTarget /
470     pm -> heapEnd - 1;
471    
472     /* Check if photon map "overflowed", i.e. exceeded its target
473     * count in the prepass; correcting the photon flux via the
474     * distribution ratio is no longer possible, as no more
475     * photons of this type will be stored, so notify the user
476     * rather than deliver incorrect results.
477     * In future we should handle this more intelligently by
478     * using the photonFlux in each photon map to individually
479     * correct the flux after distribution. */
480     if (pm -> distribRatio <= FTINY) {
481     sprintf(errmsg,
482     "%s photon map overflow in prepass, reduce -apD",
483     pmapName [t]);
484     error(INTERNAL, errmsg);
485     }
486    
487     maxDistribRatio = max(pm -> distribRatio, maxDistribRatio);
488     }
489    
490     /* Normalise distribution ratios and calculate number of photons to
491     * emit in main pass */
492     for (t = 0; t < NUM_PMAP_TYPES; t++)
493     if ((pm = photonMaps [t]))
494     pm -> distribRatio /= maxDistribRatio;
495    
496     if ((numEmit = repProgress * maxDistribRatio) < FTINY)
497     /* No photons left to distribute in main pass */
498     break;
499     }
500    
501     /* Set completion count for progress report */
502     repComplete = numEmit + repProgress;
503    
504     /* PHOTON DISTRIBUTION LOOP */
505     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
506     unsigned portCnt = 0;
507     emap.src = source + srcIdx;
508    
509     do {
510     emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
511     : NULL;
512     photonPartition [emap.src -> so -> otype] (&emap);
513    
514     if (photonRepTime) {
515     if (!passCnt)
516     sprintf(errmsg, "PREPASS %d on source %s ",
517     prePassCnt, source [srcIdx].so -> oname);
518     else
519     sprintf(errmsg, "MAIN PASS on source %s ",
520     source [srcIdx].so -> oname);
521    
522     if (emap.port) {
523     sprintf(errmsg2, "via port %s ",
524     photonPorts [portCnt].so -> oname);
525     strcat(errmsg, errmsg2);
526     }
527    
528     sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
529     strcat(errmsg, errmsg2);
530     eputs(errmsg);
531     fflush(stderr);
532     }
533    
534     for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
535     emap.partitionCnt++) {
536     double partNumEmit;
537     unsigned long partEmitCnt;
538    
539     /* Get photon origin within current source partishunn and
540     * build emission map */
541     photonOrigin [emap.src -> so -> otype] (&emap);
542     initPhotonEmission(&emap, pdfSamples);
543    
544     /* Number of photons to emit from ziss partishunn --
545     * proportional to flux; photon ray weight and scalar flux
546     * are uniform (the latter only varying in RGB). */
547     partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux;
548     partEmitCnt = (unsigned long)partNumEmit;
549    
550     /* Probabilistically account for fractional photons */
551     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
552     partEmitCnt++;
553    
554     /* Integer counter avoids FP rounding errors */
555     while (partEmitCnt--) {
556     RAY photonRay;
557    
558     /* Emit photon based on PDF and trace through scene until
559     * absorbed/leaked */
560     emitPhoton(&emap, &photonRay);
561     tracePhoton(&photonRay);
562    
563     /* Record progress */
564     repProgress++;
565    
566     if (photonRepTime > 0 &&
567     time(NULL) >= repLastTime + photonRepTime)
568     pmapDistribReport();
569 rschregle 2.3 #ifdef SIGCONT
570 greg 2.1 else signal(SIGCONT, pmapDistribReport);
571     #endif
572     }
573     }
574    
575     portCnt++;
576     } while (portCnt < numPhotonPorts);
577     }
578    
579     for (t = 0; t < NUM_PMAP_TYPES; t++)
580     if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
581     /* Double preDistrib in case a photon map is empty and redo
582     * pass 1 --> possibility of infinite loop for pathological
583     * scenes (e.g. absorbing materials) */
584     preDistrib *= 2;
585     break;
586     }
587    
588     if (t >= NUM_PMAP_TYPES) {
589     /* No empty photon maps found; now do pass 2 */
590     passCnt++;
591     if (photonRepTime)
592     eputs("\n");
593     }
594     } while (passCnt < 2);
595    
596     /* ===================================================================
597     * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
598     * =================================================================== */
599 rschregle 2.3 #ifdef SIGCONT
600     signal(SIGCONT, SIG_DFL);
601     #endif
602 greg 2.1 free(emap.samples);
603    
604     /* Set photon flux (repProgress is total num emitted) */
605     totalFlux /= repProgress;
606    
607     for (t = 0; t < NUM_PMAP_TYPES; t++)
608     if (photonMaps [t]) {
609     if (photonRepTime) {
610     sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]);
611     eputs(errmsg);
612     fflush(stderr);
613     }
614    
615     balancePhotons(photonMaps [t], &totalFlux);
616     }
617    
618     /* Precompute photon irradiance if necessary */
619     if (preCompPmap)
620     preComputeGlobal(preCompPmap);
621     }
622    
623    
624    
625     void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
626     /* Photon density estimate. Returns irradiance at ray -> rop. */
627     {
628     unsigned i;
629     PhotonSQNode *sq;
630     float r;
631     COLOR flux;
632    
633     setcolor(irrad, 0, 0, 0);
634    
635     if (!pmap -> maxGather)
636     return;
637    
638     /* Ignore sources */
639     if (ray -> ro)
640     if (islight(objptr(ray -> ro -> omod) -> otype))
641     return;
642    
643     pmap -> squeueEnd = 0;
644     findPhotons(pmap, ray);
645    
646     /* Need at least 2 photons */
647     if (pmap -> squeueEnd < 2) {
648     #ifdef PMAP_NONEFOUND
649     sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
650     ray -> ro ? ray -> ro -> oname : "<null>",
651     ray -> rop [0], ray -> rop [1], ray -> rop [2]);
652     error(WARNING, errmsg);
653     #endif
654    
655     return;
656     }
657    
658     if (pmap -> minGather == pmap -> maxGather) {
659     /* No bias compensation. Just do a plain vanilla estimate */
660     sq = pmap -> squeue + 1;
661    
662     /* Average radius between furthest two photons to improve accuracy */
663     r = max(sq -> dist, (sq + 1) -> dist);
664     r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));
665    
666     /* Skip the extra photon */
667     for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
668     getPhotonFlux(sq -> photon, flux);
669     #ifdef PMAP_EPANECHNIKOV
670     /* Apply Epanechnikov kernel to photon flux (dists are squared) */
671     scalecolor(flux, 2 * (1 - sq -> dist / r));
672     #endif
673     addcolor(irrad, flux);
674     }
675    
676     /* Divide by search area PI * r^2, 1 / PI required as ambient
677     normalisation factor */
678     scalecolor(irrad, 1 / (PI * PI * r));
679    
680     return;
681     }
682     else
683     /* Apply bias compensation to density estimate */
684     biasComp(pmap, irrad);
685     }
686    
687    
688    
689     void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
690     /* Returns precomputed photon density estimate at ray -> rop. */
691     {
692     Photon *p;
693    
694     setcolor(irrad, 0, 0, 0);
695    
696     /* Ignore sources */
697     if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
698     return;
699    
700     if ((p = find1Photon(preCompPmap, r)))
701     getPhotonFlux(p, irrad);
702     }
703    
704    
705    
706     void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
707     /* Photon volume density estimate. Returns irradiance at ray -> rop. */
708     {
709     unsigned i;
710     PhotonSQNode *sq;
711     float gecc2, r, ph;
712     COLOR flux;
713    
714     setcolor(irrad, 0, 0, 0);
715    
716     if (!pmap -> maxGather)
717     return;
718    
719     pmap -> squeueEnd = 0;
720     findPhotons(pmap, ray);
721    
722     /* Need at least 2 photons */
723     if (pmap -> squeueEnd < 2)
724     return;
725    
726     if (pmap -> minGather == pmap -> maxGather) {
727     /* No bias compensation. Just do a plain vanilla estimate */
728     gecc2 = ray -> gecc * ray -> gecc;
729     sq = pmap -> squeue + 1;
730    
731     /* Average radius between furthest two photons to improve accuracy */
732     r = max(sq -> dist, (sq + 1) -> dist);
733     r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));
734    
735     /* Skip the extra photon */
736     for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
737     /* Compute phase function for inscattering from photon */
738     if (gecc2 <= FTINY)
739     ph = 1;
740     else {
741     ph = DOT(ray -> rdir, sq -> photon -> norm) / 127;
742     ph = 1 + gecc2 - 2 * ray -> gecc * ph;
743     ph = (1 - gecc2) / (ph * sqrt(ph));
744     }
745    
746     getPhotonFlux(sq -> photon, flux);
747     scalecolor(flux, ph);
748     addcolor(irrad, flux);
749     }
750    
751     /* Divide by search volume 4 / 3 * PI * r^3 and phase function
752     normalization factor 1 / (4 * PI) */
753     scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r)));
754    
755     return;
756     }
757    
758     else
759     /* Apply bias compensation to density estimate */
760     volumeBiasComp(pmap, ray, irrad);
761     }