ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmap.c
Revision: 2.4
Committed: Thu Apr 23 20:02:04 2015 UTC (9 years ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.3: +6 -4 lines
Log Message:
Check for stale photon map skipped if octree name not initialised

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