ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmap.c
Revision: 2.7
Committed: Fri May 22 14:09:01 2015 UTC (8 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +3 -3 lines
Log Message:
Removed redundant photon map type parameter

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 greg 2.7 $Id: pmap.c,v 2.6 2015/05/21 13:54:59 greg 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 greg 2.7 static char PmapRevision [] = "$Revision: 2.6 $";
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 greg 2.7 savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv);
105 greg 2.1 }
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 rschregle 2.5 #ifdef PMAP_RUNAWAY_WARN
210 greg 2.1 error(WARNING, "runaway photon!");
211 rschregle 2.5 #endif
212 greg 2.1 return;
213     }
214 rschregle 2.5
215 greg 2.1 if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
216     return;
217    
218     if (localhit(ray, &thescene)) {
219     mod = ray -> ro -> omod;
220    
221     if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) {
222     /* Transfer ray if modifier is VOID or clipped within antimatta */
223     RAY tray;
224     photonRay(ray, &tray, PMAP_XFER, NULL);
225     tracePhoton(&tray);
226     }
227     else {
228     /* Scatter for modifier material */
229     mat = objptr(mod);
230     photonScatter [mat -> otype] (mat, ray);
231     }
232     }
233     }
234    
235    
236    
237     static void preComputeGlobal (PhotonMap *pmap)
238     /* Precompute irradiance from global photons for final gathering using
239     the first finalGather * pmap -> heapSize photons in the heap. Returns
240     new heap with precomputed photons. */
241     {
242     unsigned long i, nuHeapSize;
243     unsigned j;
244     Photon *nuHeap, *p;
245     COLOR irrad;
246     RAY ray;
247     float nuMinPos [3], nuMaxPos [3];
248    
249     repComplete = nuHeapSize = finalGather * pmap -> heapSize;
250    
251     if (photonRepTime) {
252     sprintf(errmsg,
253     "Precomputing irradiance for %ld global photons...\n",
254     nuHeapSize);
255     eputs(errmsg);
256     fflush(stderr);
257     }
258    
259     p = nuHeap = (Photon*)malloc(nuHeapSize * sizeof(Photon));
260     if (!nuHeap)
261     error(USER, "can't allocate photon heap");
262    
263     for (j = 0; j <= 2; j++) {
264     nuMinPos [j] = FHUGE;
265     nuMaxPos [j] = -FHUGE;
266     }
267    
268     /* Record start time, baby */
269     repStartTime = time(NULL);
270 rschregle 2.3 #ifdef SIGCONT
271     signal(SIGCONT, pmapPreCompReport);
272     #endif
273 greg 2.1 repProgress = 0;
274 greg 2.2 memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon));
275 greg 2.1
276     for (i = 0, p = nuHeap; i < nuHeapSize; i++, p++) {
277     ray.ro = NULL;
278     VCOPY(ray.rop, p -> pos);
279    
280     /* Update min and max positions & set ray normal */
281     for (j = 0; j < 3; j++) {
282     if (p -> pos [j] < nuMinPos [j]) nuMinPos [j] = p -> pos [j];
283     if (p -> pos [j] > nuMaxPos [j]) nuMaxPos [j] = p -> pos [j];
284     ray.ron [j] = p -> norm [j] / 127.0;
285     }
286    
287     photonDensity(pmap, &ray, irrad);
288     setPhotonFlux(p, irrad);
289     repProgress++;
290    
291     if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
292     pmapPreCompReport();
293 rschregle 2.3 #ifdef SIGCONT
294     else signal(SIGCONT, pmapPreCompReport);
295     #endif
296 greg 2.1 }
297    
298 rschregle 2.3 #ifdef SIGCONT
299     signal(SIGCONT, SIG_DFL);
300     #endif
301 greg 2.1
302     /* Replace & rebuild heap */
303     free(pmap -> heap);
304     pmap -> heap = nuHeap;
305     pmap -> heapSize = pmap -> heapEnd = nuHeapSize;
306     VCOPY(pmap -> minPos, nuMinPos);
307     VCOPY(pmap -> maxPos, nuMaxPos);
308    
309     if (photonRepTime) {
310     eputs("Rebuilding global photon heap...\n");
311     fflush(stderr);
312     }
313    
314     balancePhotons(pmap, NULL);
315     }
316    
317    
318    
319     void distribPhotons (PhotonMap **pmaps)
320     {
321     EmissionMap emap;
322     char errmsg2 [128];
323     unsigned t, srcIdx, passCnt = 0, prePassCnt = 0;
324     double totalFlux = 0;
325     PhotonMap *pm;
326    
327     for (t = 0; t < NUM_PMAP_TYPES && !photonMaps [t]; t++);
328     if (t >= NUM_PMAP_TYPES)
329     error(USER, "no photon maps defined");
330    
331     if (!nsources)
332     error(USER, "no light sources");
333    
334     /* ===================================================================
335     * INITIALISATION - Set up emission and scattering funcs
336     * =================================================================== */
337     emap.samples = NULL;
338     emap.maxPartitions = MAXSPART;
339     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
340     if (!emap.partitions)
341     error(INTERNAL, "can't allocate source partitions");
342    
343     /* Initialise all defined photon maps */
344     for (t = 0; t < NUM_PMAP_TYPES; t++)
345     initPhotonMap(photonMaps [t], t);
346    
347     initPhotonEmissionFuncs();
348     initPhotonScatterFuncs();
349    
350     /* Get photon ports if specified */
351     if (ambincl == 1)
352     getPhotonPorts();
353    
354     /* Get photon sensor modifiers */
355     getPhotonSensors(photonSensorList);
356    
357     /* Seed RNGs for photon distribution */
358     pmapSeed(randSeed, partState);
359     pmapSeed(randSeed, emitState);
360     pmapSeed(randSeed, cntState);
361     pmapSeed(randSeed, mediumState);
362     pmapSeed(randSeed, scatterState);
363     pmapSeed(randSeed, rouletteState);
364    
365     if (photonRepTime)
366     eputs("\n");
367    
368     /* ===================================================================
369     * FLUX INTEGRATION - Get total photon flux from light sources
370     * =================================================================== */
371     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
372     unsigned portCnt = 0;
373     emap.src = source + srcIdx;
374    
375     do {
376     emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
377     : NULL;
378     photonPartition [emap.src -> so -> otype] (&emap);
379    
380     if (photonRepTime) {
381     sprintf(errmsg, "Integrating flux from source %s ",
382     source [srcIdx].so -> oname);
383    
384     if (emap.port) {
385     sprintf(errmsg2, "via port %s ",
386     photonPorts [portCnt].so -> oname);
387     strcat(errmsg, errmsg2);
388     }
389    
390     sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
391     strcat(errmsg, errmsg2);
392     eputs(errmsg);
393     fflush(stderr);
394     }
395    
396     for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
397     emap.partitionCnt++) {
398     initPhotonEmission(&emap, pdfSamples);
399     totalFlux += colorAvg(emap.partFlux);
400     }
401    
402     portCnt++;
403     } while (portCnt < numPhotonPorts);
404     }
405    
406     if (totalFlux < FTINY)
407     error(USER, "zero flux from light sources");
408    
409     /* Record start time and enable progress report signal handler */
410     repStartTime = time(NULL);
411 rschregle 2.3 #ifdef SIGCONT
412     signal(SIGCONT, pmapDistribReport);
413     #endif
414 greg 2.1 repProgress = prePassCnt = 0;
415    
416     if (photonRepTime)
417     eputs("\n");
418    
419     /* ===================================================================
420     * 2-PASS PHOTON DISTRIBUTION
421     * Pass 1 (pre): emit fraction of target photon count
422     * Pass 2 (main): based on outcome of pass 1, estimate remaining number
423     * of photons to emit to approximate target count
424     * =================================================================== */
425     do {
426     double numEmit;
427    
428     if (!passCnt) {
429     /* INIT PASS 1 */
430     /* Skip if no photons contributed after sufficient iterations; make
431     * it clear to user which photon maps are missing so (s)he can
432     * check the scene geometry and materials */
433     if (++prePassCnt > maxPreDistrib) {
434     sprintf(errmsg, "too many prepasses");
435    
436     for (t = 0; t < NUM_PMAP_TYPES; t++)
437     if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
438     sprintf(errmsg2, ", no %s photons stored", pmapName [t]);
439     strcat(errmsg, errmsg2);
440     }
441    
442     error(USER, errmsg);
443     break;
444     }
445    
446     /* Num to emit is fraction of minimum target count */
447     numEmit = FHUGE;
448    
449     for (t = 0; t < NUM_PMAP_TYPES; t++)
450     if (photonMaps [t])
451     numEmit = min(photonMaps [t] -> distribTarget, numEmit);
452    
453     numEmit *= preDistrib;
454     }
455    
456     else {
457     /* INIT PASS 2 */
458     /* Based on the outcome of the predistribution we can now estimate
459     * how many more photons we have to emit for each photon map to
460     * meet its respective target count. This value is clamped to 0 in
461     * case the target has already been exceeded in the pass 1. Note
462     * repProgress is the number of photons emitted thus far, while
463     * heapEnd is the number of photons stored in each photon map. */
464     double maxDistribRatio = 0;
465    
466     /* Set the distribution ratio for each map; this indicates how many
467     * photons of each respective type are stored per emitted photon,
468     * and is used as probability for storing a photon by addPhoton().
469     * Since this biases the photon density, addPhoton() promotes the
470     * flux of stored photons to compensate. */
471     for (t = 0; t < NUM_PMAP_TYPES; t++)
472     if ((pm = photonMaps [t])) {
473     pm -> distribRatio = (double)pm -> distribTarget /
474     pm -> heapEnd - 1;
475    
476     /* Check if photon map "overflowed", i.e. exceeded its target
477     * count in the prepass; correcting the photon flux via the
478     * distribution ratio is no longer possible, as no more
479     * photons of this type will be stored, so notify the user
480     * rather than deliver incorrect results.
481     * In future we should handle this more intelligently by
482     * using the photonFlux in each photon map to individually
483     * correct the flux after distribution. */
484     if (pm -> distribRatio <= FTINY) {
485     sprintf(errmsg,
486     "%s photon map overflow in prepass, reduce -apD",
487     pmapName [t]);
488     error(INTERNAL, errmsg);
489     }
490    
491     maxDistribRatio = max(pm -> distribRatio, maxDistribRatio);
492     }
493    
494     /* Normalise distribution ratios and calculate number of photons to
495     * emit in main pass */
496     for (t = 0; t < NUM_PMAP_TYPES; t++)
497     if ((pm = photonMaps [t]))
498     pm -> distribRatio /= maxDistribRatio;
499    
500     if ((numEmit = repProgress * maxDistribRatio) < FTINY)
501     /* No photons left to distribute in main pass */
502     break;
503     }
504    
505     /* Set completion count for progress report */
506     repComplete = numEmit + repProgress;
507    
508     /* PHOTON DISTRIBUTION LOOP */
509     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
510     unsigned portCnt = 0;
511     emap.src = source + srcIdx;
512 greg 2.6
513 greg 2.1 do {
514     emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
515     : NULL;
516     photonPartition [emap.src -> so -> otype] (&emap);
517    
518     if (photonRepTime) {
519     if (!passCnt)
520     sprintf(errmsg, "PREPASS %d on source %s ",
521     prePassCnt, source [srcIdx].so -> oname);
522     else
523     sprintf(errmsg, "MAIN PASS on source %s ",
524     source [srcIdx].so -> oname);
525    
526     if (emap.port) {
527     sprintf(errmsg2, "via port %s ",
528     photonPorts [portCnt].so -> oname);
529     strcat(errmsg, errmsg2);
530     }
531    
532     sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
533     strcat(errmsg, errmsg2);
534     eputs(errmsg);
535     fflush(stderr);
536     }
537    
538     for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
539     emap.partitionCnt++) {
540     double partNumEmit;
541     unsigned long partEmitCnt;
542    
543     /* Get photon origin within current source partishunn and
544     * build emission map */
545     photonOrigin [emap.src -> so -> otype] (&emap);
546     initPhotonEmission(&emap, pdfSamples);
547    
548     /* Number of photons to emit from ziss partishunn --
549     * proportional to flux; photon ray weight and scalar flux
550     * are uniform (the latter only varying in RGB). */
551     partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux;
552     partEmitCnt = (unsigned long)partNumEmit;
553    
554     /* Probabilistically account for fractional photons */
555     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
556     partEmitCnt++;
557    
558     /* Integer counter avoids FP rounding errors */
559     while (partEmitCnt--) {
560     RAY photonRay;
561    
562     /* Emit photon based on PDF and trace through scene until
563     * absorbed/leaked */
564     emitPhoton(&emap, &photonRay);
565     tracePhoton(&photonRay);
566    
567     /* Record progress */
568     repProgress++;
569    
570     if (photonRepTime > 0 &&
571     time(NULL) >= repLastTime + photonRepTime)
572     pmapDistribReport();
573 rschregle 2.3 #ifdef SIGCONT
574 greg 2.1 else signal(SIGCONT, pmapDistribReport);
575     #endif
576     }
577     }
578    
579     portCnt++;
580     } while (portCnt < numPhotonPorts);
581     }
582    
583     for (t = 0; t < NUM_PMAP_TYPES; t++)
584     if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
585     /* Double preDistrib in case a photon map is empty and redo
586     * pass 1 --> possibility of infinite loop for pathological
587     * scenes (e.g. absorbing materials) */
588     preDistrib *= 2;
589     break;
590     }
591    
592     if (t >= NUM_PMAP_TYPES) {
593     /* No empty photon maps found; now do pass 2 */
594     passCnt++;
595     if (photonRepTime)
596     eputs("\n");
597     }
598     } while (passCnt < 2);
599    
600     /* ===================================================================
601     * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
602     * =================================================================== */
603 rschregle 2.3 #ifdef SIGCONT
604     signal(SIGCONT, SIG_DFL);
605     #endif
606 greg 2.1 free(emap.samples);
607    
608     /* Set photon flux (repProgress is total num emitted) */
609     totalFlux /= repProgress;
610    
611     for (t = 0; t < NUM_PMAP_TYPES; t++)
612     if (photonMaps [t]) {
613     if (photonRepTime) {
614     sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]);
615     eputs(errmsg);
616     fflush(stderr);
617     }
618    
619     balancePhotons(photonMaps [t], &totalFlux);
620     }
621    
622     /* Precompute photon irradiance if necessary */
623     if (preCompPmap)
624     preComputeGlobal(preCompPmap);
625     }
626    
627    
628    
629     void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
630     /* Photon density estimate. Returns irradiance at ray -> rop. */
631     {
632     unsigned i;
633     PhotonSQNode *sq;
634     float r;
635     COLOR flux;
636    
637     setcolor(irrad, 0, 0, 0);
638    
639     if (!pmap -> maxGather)
640     return;
641    
642     /* Ignore sources */
643     if (ray -> ro)
644     if (islight(objptr(ray -> ro -> omod) -> otype))
645     return;
646    
647     pmap -> squeueEnd = 0;
648     findPhotons(pmap, ray);
649    
650     /* Need at least 2 photons */
651     if (pmap -> squeueEnd < 2) {
652     #ifdef PMAP_NONEFOUND
653     sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
654     ray -> ro ? ray -> ro -> oname : "<null>",
655     ray -> rop [0], ray -> rop [1], ray -> rop [2]);
656     error(WARNING, errmsg);
657     #endif
658    
659     return;
660     }
661    
662     if (pmap -> minGather == pmap -> maxGather) {
663     /* No bias compensation. Just do a plain vanilla estimate */
664     sq = pmap -> squeue + 1;
665    
666     /* Average radius between furthest two photons to improve accuracy */
667     r = max(sq -> dist, (sq + 1) -> dist);
668     r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));
669    
670     /* Skip the extra photon */
671     for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
672     getPhotonFlux(sq -> photon, flux);
673     #ifdef PMAP_EPANECHNIKOV
674     /* Apply Epanechnikov kernel to photon flux (dists are squared) */
675     scalecolor(flux, 2 * (1 - sq -> dist / r));
676     #endif
677     addcolor(irrad, flux);
678     }
679    
680     /* Divide by search area PI * r^2, 1 / PI required as ambient
681     normalisation factor */
682     scalecolor(irrad, 1 / (PI * PI * r));
683    
684     return;
685     }
686     else
687     /* Apply bias compensation to density estimate */
688     biasComp(pmap, irrad);
689     }
690    
691    
692    
693     void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
694     /* Returns precomputed photon density estimate at ray -> rop. */
695     {
696     Photon *p;
697    
698     setcolor(irrad, 0, 0, 0);
699    
700     /* Ignore sources */
701     if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
702     return;
703    
704     if ((p = find1Photon(preCompPmap, r)))
705     getPhotonFlux(p, irrad);
706     }
707    
708    
709    
710     void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
711     /* Photon volume density estimate. Returns irradiance at ray -> rop. */
712     {
713     unsigned i;
714     PhotonSQNode *sq;
715     float gecc2, r, ph;
716     COLOR flux;
717    
718     setcolor(irrad, 0, 0, 0);
719    
720     if (!pmap -> maxGather)
721     return;
722    
723     pmap -> squeueEnd = 0;
724     findPhotons(pmap, ray);
725    
726     /* Need at least 2 photons */
727     if (pmap -> squeueEnd < 2)
728     return;
729    
730     if (pmap -> minGather == pmap -> maxGather) {
731     /* No bias compensation. Just do a plain vanilla estimate */
732     gecc2 = ray -> gecc * ray -> gecc;
733     sq = pmap -> squeue + 1;
734    
735     /* Average radius between furthest two photons to improve accuracy */
736     r = max(sq -> dist, (sq + 1) -> dist);
737     r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));
738    
739     /* Skip the extra photon */
740     for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
741     /* Compute phase function for inscattering from photon */
742     if (gecc2 <= FTINY)
743     ph = 1;
744     else {
745     ph = DOT(ray -> rdir, sq -> photon -> norm) / 127;
746     ph = 1 + gecc2 - 2 * ray -> gecc * ph;
747     ph = (1 - gecc2) / (ph * sqrt(ph));
748     }
749    
750     getPhotonFlux(sq -> photon, flux);
751     scalecolor(flux, ph);
752     addcolor(irrad, flux);
753     }
754    
755     /* Divide by search volume 4 / 3 * PI * r^3 and phase function
756     normalization factor 1 / (4 * PI) */
757     scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r)));
758    
759     return;
760     }
761    
762     else
763     /* Apply bias compensation to density estimate */
764     volumeBiasComp(pmap, ray, irrad);
765     }