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

# Content
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 (c) Lucerne University of Applied Sciences and Arts,
8 supported by the Swiss National Science Foundation (SNSF, #147053)
9 ==================================================================
10
11 $Id: pmap.c,v 2.6 2015/05/21 13:54:59 greg Exp $
12 */
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 static char PmapRevision [] = "$Revision: 2.6 $";
32
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 if (pm -> fileName && octname &&
69 !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) &&
70 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, 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 #ifdef PMAP_RUNAWAY_WARN
210 error(WARNING, "runaway photon!");
211 #endif
212 return;
213 }
214
215 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 #ifdef SIGCONT
271 signal(SIGCONT, pmapPreCompReport);
272 #endif
273 repProgress = 0;
274 memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon));
275
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 #ifdef SIGCONT
294 else signal(SIGCONT, pmapPreCompReport);
295 #endif
296 }
297
298 #ifdef SIGCONT
299 signal(SIGCONT, SIG_DFL);
300 #endif
301
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 #ifdef SIGCONT
412 signal(SIGCONT, pmapDistribReport);
413 #endif
414 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
513 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 #ifdef SIGCONT
574 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 #ifdef SIGCONT
604 signal(SIGCONT, SIG_DFL);
605 #endif
606 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 }