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 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

# 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 Lucerne University of Applied Sciences & Arts
8 ==================================================================
9
10 $Id: pmap.c,v 2.2 2015/04/21 19:16:51 greg Exp $
11 */
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 static char PmapRevision [] = "$Revision: 2.2 $";
31
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 #ifdef SIGCONT
267 signal(SIGCONT, pmapPreCompReport);
268 #endif
269 repProgress = 0;
270 memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon));
271
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 #ifdef SIGCONT
290 else signal(SIGCONT, pmapPreCompReport);
291 #endif
292 }
293
294 #ifdef SIGCONT
295 signal(SIGCONT, SIG_DFL);
296 #endif
297
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 #ifdef SIGCONT
408 signal(SIGCONT, pmapDistribReport);
409 #endif
410 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 #ifdef SIGCONT
570 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 #ifdef SIGCONT
600 signal(SIGCONT, SIG_DFL);
601 #endif
602 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 }