ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapdata.c
Revision: 2.2
Committed: Wed Apr 22 15:44:57 2015 UTC (9 years ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.1: +3 -1 lines
Log Message:
Conditionally disabled warning about clamped max photon lookup radius.

File Contents

# User Rev Content
1 greg 2.1 /*
2     ==================================================================
3     Photon map data structures and kd-tree handling
4    
5     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
6     (c) Fraunhofer Institute for Solar Energy Systems,
7     Lucerne University of Applied Sciences & Arts
8     ==================================================================
9    
10 rschregle 2.2 $Id: pmapdata.c,v 4.31 2015/04/22 15:37:07 taschreg Exp taschreg $
11 greg 2.1 */
12    
13    
14    
15     #include "pmap.h"
16     #include "pmaprand.h"
17     #include "pmapmat.h"
18     #include "otypes.h"
19     #include "source.h"
20     #include "rcontrib.h"
21    
22    
23    
24     PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
25     NULL, NULL, NULL, NULL, NULL, NULL
26     };
27    
28    
29    
30     void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
31     /* Init photon map 'n' stuff... */
32     {
33     if (!pmap)
34     return;
35    
36     pmap -> heapSize = pmap -> heapEnd = 0;
37     pmap -> heap = NULL;
38     pmap -> squeue = NULL;
39     pmap -> biasCompHist = NULL;
40     pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE;
41     pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE;
42     pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0;
43     pmap -> gatherTolerance = gatherTolerance;
44     pmap -> minError = pmap -> maxError = pmap -> rmsError = 0;
45     pmap -> numDensity = 0;
46     pmap -> distribRatio = 1;
47     pmap -> type = t;
48    
49     /* Init local RNG state */
50     pmap -> randState [0] = 10243;
51     pmap -> randState [1] = 39829;
52     pmap -> randState [2] = 9433;
53     /* pmapSeed(25999, pmap -> randState); */
54     pmapSeed(randSeed, pmap -> randState);
55    
56     /* Set up type-specific photon lookup callback */
57     pmap -> lookup = pmapLookup [t];
58    
59     pmap -> primary = NULL;
60     pmap -> primarySize = pmap -> primaryEnd = 0;
61     }
62    
63    
64    
65     const PhotonPrimary* addPhotonPrimary (PhotonMap *pmap, const RAY *ray)
66     {
67     PhotonPrimary *prim = NULL;
68    
69     if (!pmap || !ray)
70     return NULL;
71    
72     /* Check if last primary ray has spawned photons (srcIdx >= 0, see
73     * addPhoton()), in which case we keep it and allocate a new one;
74     * otherwise we overwrite the unused entry */
75     if (pmap -> primary && pmap -> primary [pmap -> primaryEnd].srcIdx >= 0)
76     pmap -> primaryEnd++;
77    
78     if (!pmap -> primarySize || pmap -> primaryEnd >= pmap -> primarySize) {
79     /* Allocate/enlarge array */
80     pmap -> primarySize += pmap -> heapSizeInc;
81    
82     /* Counter wraparound? */
83     if (pmap -> primarySize < pmap -> heapSizeInc)
84     error(INTERNAL, "photon primary overflow");
85    
86     pmap -> primary = (PhotonPrimary *)realloc(pmap -> primary,
87     sizeof(PhotonPrimary) *
88     pmap -> primarySize);
89     if (!pmap -> primary)
90     error(USER, "can't allocate photon primaries");
91     }
92    
93     prim = pmap -> primary + pmap -> primaryEnd;
94    
95     /* Mark unused with negative source index until path spawns a photon (see
96     * addPhoton()) */
97     prim -> srcIdx = -1;
98    
99     /* Reverse incident direction to point to light source */
100     prim -> dir [0] = -ray -> rdir [0];
101     prim -> dir [1] = -ray -> rdir [1];
102     prim -> dir [2] = -ray -> rdir [2];
103    
104     VCOPY(prim -> org, ray -> rorg);
105    
106     return prim;
107     }
108    
109    
110    
111     const Photon* addPhoton (PhotonMap* pmap, const RAY* ray)
112     {
113     unsigned i;
114     Photon* photon = NULL;
115     COLOR photonFlux;
116    
117     /* Account for distribution ratio */
118     if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio)
119     return NULL;
120    
121     /* Don't store on sources */
122     if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
123     return NULL;
124    
125     #if 0
126     if (inContribPmap(pmap)) {
127     /* Adding contribution photon */
128     if (ray -> parent && ray -> parent -> rtype & PRIMARY)
129     /* Add primary photon ray if parent is primary; by putting this
130     * here and checking the ray's immediate parent, we only add
131     * primaries that actually contribute photons, and we only add them
132     * once. */
133     addPhotonPrimary(pmap, ray -> parent);
134    
135     /* Save index to primary ray (remains unchanged if primary already in
136     * array) */
137     primary = pmap -> primaryEnd;
138     }
139     #endif
140    
141     #ifdef PMAP_ROI
142     /* Store photon if within region of interest -- for egg-spurtz only! */
143     if (ray -> rop [0] >= pmapROI [0] && ray -> rop [0] <= pmapROI [1] &&
144     ray -> rop [1] >= pmapROI [2] && ray -> rop [1] <= pmapROI [3] &&
145     ray -> rop [2] >= pmapROI [4] && ray -> rop [2] <= pmapROI [5])
146     #endif
147     {
148     if (pmap -> heapEnd >= pmap -> heapSize) {
149     /* Enlarge heap */
150     pmap -> heapSize += pmap -> heapSizeInc;
151    
152     /* Counter wraparound? */
153     if (pmap -> heapSize < pmap -> heapSizeInc)
154     error(INTERNAL, "photon heap overflow");
155    
156     pmap -> heap = (Photon *)realloc(pmap -> heap,
157     sizeof(Photon) * pmap -> heapSize);
158     if (!pmap -> heap)
159     error(USER, "can't allocate photon heap");
160     }
161    
162     photon = pmap -> heap + pmap -> heapEnd++;
163    
164     /* Adjust flux according to distribution ratio and ray weight */
165     copycolor(photonFlux, ray -> rcol);
166     scalecolor(photonFlux,
167     ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
168     : 1));
169     setPhotonFlux(photon, photonFlux);
170    
171     /* Set photon position and flags */
172     VCOPY(photon -> pos, ray -> rop);
173     photon -> flags = PMAP_CAUSTICRAY(ray) ? PMAP_CAUST_BIT : 0;
174    
175     /* Set primary ray index and mark as used for contrib photons */
176     if (isContribPmap(pmap)) {
177     photon -> primary = pmap -> primaryEnd;
178     pmap -> primary [pmap -> primaryEnd].srcIdx = ray -> rsrc;
179     }
180     else photon -> primary = 0;
181    
182     /* Update min and max positions & set normal */
183     for (i = 0; i <= 2; i++) {
184     if (photon -> pos [i] < pmap -> minPos [i])
185     pmap -> minPos [i] = photon -> pos [i];
186     if (photon -> pos [i] > pmap -> maxPos [i])
187     pmap -> maxPos [i] = photon -> pos [i];
188     photon -> norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
189     : ray -> ron [i]);
190     }
191     }
192    
193     return photon;
194     }
195    
196    
197    
198     static void nearestNeighbours (PhotonMap* pmap, const float pos [3],
199     const float norm [3], unsigned long node)
200     /* Recursive part of findPhotons(..).
201     Note that all heap and priority queue index handling is 1-based, but
202     accesses to the arrays are 0-based! */
203     {
204     Photon* p = &pmap -> heap [node - 1];
205     unsigned i, j;
206     /* Signed distance to current photon's splitting plane */
207     float d = pos [photonDiscr(*p)] - p -> pos [photonDiscr(*p)],
208     d2 = d * d;
209     PhotonSQNode* sq = pmap -> squeue;
210     const unsigned sqSize = pmap -> squeueSize;
211     float dv [3];
212    
213     /* Search subtree closer to pos first; exclude other subtree if the
214     distance to the splitting plane is greater than maxDist */
215     if (d < 0) {
216     if (node << 1 <= pmap -> heapSize)
217     nearestNeighbours(pmap, pos, norm, node << 1);
218     if (d2 < pmap -> maxDist && node << 1 < pmap -> heapSize)
219     nearestNeighbours(pmap, pos, norm, (node << 1) + 1);
220     }
221     else {
222     if (node << 1 < pmap -> heapSize)
223     nearestNeighbours(pmap, pos, norm, (node << 1) + 1);
224     if (d2 < pmap -> maxDist && node << 1 <= pmap -> heapSize)
225     nearestNeighbours(pmap, pos, norm, node << 1);
226     }
227    
228     /* Reject photon if normal faces away (ignored for volume photons) */
229     if (norm && DOT(norm, p -> norm) <= 0)
230     return;
231    
232     if (isContribPmap(pmap) && pmap -> srcContrib) {
233     /* Lookup in contribution photon map */
234     OBJREC *srcMod;
235     const int srcIdx = photonSrcIdx(pmap, p);
236    
237     if (srcIdx < 0 || srcIdx >= nsources)
238     error(INTERNAL, "invalid light source index in photon map");
239    
240     srcMod = objptr(source [srcIdx].so -> omod);
241    
242     /* Reject photon if contributions from light source which emitted it
243     * are not sought */
244     if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
245     return;
246    
247     /* Reject non-caustic photon if lookup for caustic contribs */
248     if (pmap -> lookupFlags & PMAP_CAUST_BIT & ~p -> flags)
249     return;
250     }
251    
252     /* Squared distance to current photon */
253     dv [0] = pos [0] - p -> pos [0];
254     dv [1] = pos [1] - p -> pos [1];
255     dv [2] = pos [2] - p -> pos [2];
256     d2 = DOT(dv, dv);
257    
258     /* Accept photon if closer than current max dist & add to priority queue */
259     if (d2 < pmap -> maxDist) {
260     if (pmap -> squeueEnd < sqSize) {
261     /* Priority queue not full; append photon and restore heap */
262     i = ++pmap -> squeueEnd;
263    
264     while (i > 1 && sq [(i >> 1) - 1].dist <= d2) {
265     sq [i - 1].photon = sq [(i >> 1) - 1].photon;
266     sq [i - 1].dist = sq [(i >> 1) - 1].dist;
267     i >>= 1;
268     }
269    
270     sq [--i].photon = p;
271     sq [i].dist = d2;
272     /* Update maxDist if we've just filled the queue */
273     if (pmap -> squeueEnd >= pmap -> squeueSize)
274     pmap -> maxDist = sq [0].dist;
275     }
276     else {
277     /* Priority queue full; replace maximum, restore heap, and
278     update maxDist */
279     i = 1;
280    
281     while (i <= sqSize >> 1) {
282     j = i << 1;
283     if (j < sqSize && sq [j - 1].dist < sq [j].dist)
284     j++;
285     if (d2 >= sq [j - 1].dist)
286     break;
287     sq [i - 1].photon = sq [j - 1].photon;
288     sq [i - 1].dist = sq [j - 1].dist;
289     i = j;
290     }
291    
292     sq [--i].photon = p;
293     sq [i].dist = d2;
294     pmap -> maxDist = sq [0].dist;
295     }
296     }
297     }
298    
299    
300    
301     /* Dynamic max photon search radius increase and reduction factors */
302     #define PMAP_MAXDIST_INC 4
303     #define PMAP_MAXDIST_DEC 0.9
304    
305     /* Num successful lookups before reducing in max search radius */
306     #define PMAP_MAXDIST_CNT 1000
307    
308     /* Threshold below which we assume increasing max radius won't help */
309     #define PMAP_SHORT_LOOKUP_THRESH 1
310    
311     void findPhotons (PhotonMap* pmap, const RAY* ray)
312     {
313     float pos [3], norm [3];
314     int redo = 0;
315    
316     if (!pmap -> squeue) {
317     /* Lazy init priority queue */
318     pmap -> squeueSize = pmap -> maxGather + 1;
319     pmap -> squeue = (PhotonSQNode*)malloc(pmap -> squeueSize *
320     sizeof(PhotonSQNode));
321     if (!pmap -> squeue)
322     error(USER, "can't allocate photon priority queue");
323    
324     pmap -> minGathered = pmap -> maxGather;
325     pmap -> maxGathered = pmap -> minGather;
326     pmap -> totalGathered = 0;
327     pmap -> numLookups = pmap -> numShortLookups = 0;
328     pmap -> shortLookupPct = 0;
329     pmap -> minError = FHUGE;
330     pmap -> maxError = -FHUGE;
331     pmap -> rmsError = 0;
332     /* Maximum search radius limit based on avg photon distance to
333     * centre of gravity */
334     pmap -> maxDist0 = pmap -> maxDistLimit =
335     maxDistCoeff * pmap -> squeueSize * pmap -> CoGdist /
336     pmap -> heapSize;
337     }
338    
339     do {
340     pmap -> squeueEnd = 0;
341     pmap -> maxDist = pmap -> maxDist0;
342    
343     /* Search position is ray -> rorg for volume photons, since we have no
344     intersection point. Normals are ignored -- these are incident
345     directions). */
346     if (isVolumePmap(pmap)) {
347     VCOPY(pos, ray -> rorg);
348     nearestNeighbours(pmap, pos, NULL, 1);
349     }
350     else {
351     VCOPY(pos, ray -> rop);
352     VCOPY(norm, ray -> ron);
353     nearestNeighbours(pmap, pos, norm, 1);
354     }
355    
356     if (pmap -> squeueEnd < pmap -> squeueSize * pmap -> gatherTolerance) {
357     /* Short lookup; too few photons found */
358     if (pmap -> squeueEnd > PMAP_SHORT_LOOKUP_THRESH) {
359     /* Ignore short lookups which return fewer than
360     * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there
361     * really are no photons in the vicinity, and increasing the max
362     * search radius therefore won't help */
363     #ifdef PMAP_LOOKUP_WARN
364     sprintf(errmsg,
365     "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
366     pmap -> squeueEnd, pmap -> squeueSize,
367     pmapName [pmap -> type], pos [0], pos [1], pos [2],
368     ray -> ro ? ray -> ro -> oname : "<null>");
369     error(WARNING, errmsg);
370     #endif
371    
372     if (pmap -> maxDist0 < pmap -> maxDistLimit) {
373     /* Increase max search radius if below limit & redo search */
374     pmap -> maxDist0 *= PMAP_MAXDIST_INC;
375     #ifdef PMAP_LOOKUP_REDO
376     redo = 1;
377     #endif
378    
379     #ifdef PMAP_LOOKUP_WARN
380     sprintf(errmsg,
381     redo ? "restarting photon lookup with max radius %.1e"
382     : "max photon lookup radius adjusted to %.1e",
383     pmap -> maxDist0);
384     error(WARNING, errmsg);
385     #endif
386     }
387 rschregle 2.2 #ifdef PMAP_LOOKUP_REDO
388 greg 2.1 else {
389     sprintf(errmsg, "max photon lookup radius clamped to %.1e",
390     pmap -> maxDist0);
391     error(WARNING, errmsg);
392     }
393 rschregle 2.2 #endif
394 greg 2.1 }
395    
396     /* Reset successful lookup counter */
397     pmap -> numLookups = 0;
398     }
399     else {
400     /* Increment successful lookup counter and reduce max search radius if
401     * wraparound */
402     pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT;
403     if (!pmap -> numLookups)
404     pmap -> maxDist0 *= PMAP_MAXDIST_DEC;
405    
406     redo = 0;
407     }
408     } while (redo);
409     }
410    
411    
412    
413     static void nearest1Neighbour (PhotonMap *pmap, const float pos [3],
414     const float norm [3], Photon **photon,
415     unsigned long node)
416     /* Recursive part of find1Photon(..).
417     Note that all heap index handling is 1-based, but accesses to the
418     arrays are 0-based! */
419     {
420     Photon *p = pmap -> heap + node - 1;
421     /* Signed distance to current photon's splitting plane */
422     float d = pos [photonDiscr(*p)] - p -> pos [photonDiscr(*p)],
423     d2 = d * d;
424     float dv [3];
425    
426     /* Search subtree closer to pos first; exclude other subtree if the
427     distance to the splitting plane is greater than maxDist */
428     if (d < 0) {
429     if (node << 1 <= pmap -> heapSize)
430     nearest1Neighbour(pmap, pos, norm, photon, node << 1);
431     if (d2 < pmap -> maxDist && node << 1 < pmap -> heapSize)
432     nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
433     }
434     else {
435     if (node << 1 < pmap -> heapSize)
436     nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
437     if (d2 < pmap -> maxDist && node << 1 <= pmap -> heapSize)
438     nearest1Neighbour(pmap, pos, norm, photon, node << 1);
439     }
440    
441     /* Squared distance to current photon */
442     dv [0] = pos [0] - p -> pos [0];
443     dv [1] = pos [1] - p -> pos [1];
444     dv [2] = pos [2] - p -> pos [2];
445     d2 = DOT(dv, dv);
446    
447     if (d2 < pmap -> maxDist && DOT(norm, p -> norm) > 0) {
448     /* Closest photon so far with similar normal */
449     pmap -> maxDist = d2;
450     *photon = p;
451     }
452     }
453    
454    
455    
456     Photon* find1Photon (PhotonMap *pmap, const RAY* ray)
457     {
458     float fpos [3], norm [3];
459     Photon* photon = NULL;
460    
461     VCOPY(fpos, ray -> rop);
462     VCOPY(norm, ray -> ron);
463     pmap -> maxDist = thescene.cusize;
464     nearest1Neighbour(pmap, fpos, norm, &photon, 1);
465    
466     return photon;
467     }
468    
469    
470    
471     static unsigned long medianPartition (const Photon* heap,
472     unsigned long* heapIdx,
473     unsigned long* heapXdi,
474     unsigned long left,
475     unsigned long right, unsigned dim)
476     /* Returns index to median in heap from indices left to right
477     (inclusive) in dimension dim. The heap is partitioned relative to
478     median using a quicksort algorithm. The heap indices in heapIdx are
479     sorted rather than the heap itself. */
480     {
481     register const float* p;
482     const unsigned long n = right - left + 1;
483     register unsigned long l, r, lg2, n2, m;
484     register unsigned d;
485    
486     /* Round down n to nearest power of 2 */
487     for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2);
488     n2 = 1 << lg2;
489    
490     /* Determine median position; this takes into account the fact that
491     only the last level in the heap can be partially empty, and that
492     it fills from left to right */
493     m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1));
494    
495     while (right > left) {
496     /* Pivot node */
497     p = heap [heapIdx [right]].pos;
498     l = left;
499     r = right - 1;
500    
501     /* l & r converge, swapping elements out of order with respect to
502     pivot node. Identical keys are resolved by cycling through
503     dim. The convergence point is then the pivot's position. */
504     do {
505     while (l <= r) {
506     d = dim;
507    
508     while (heap [heapIdx [l]].pos [d] == p [d]) {
509     d = (d + 1) % 3;
510    
511     if (d == dim) {
512     /* Ignore dupes? */
513     error(WARNING, "duplicate keys in photon heap");
514     l++;
515     break;
516     }
517     }
518    
519     if (heap [heapIdx [l]].pos [d] < p [d])
520     l++;
521     else break;
522     }
523    
524     while (r > l) {
525     d = dim;
526    
527     while (heap [heapIdx [r]].pos [d] == p [d]) {
528     d = (d + 1) % 3;
529    
530     if (d == dim) {
531     /* Ignore dupes? */
532     error(WARNING, "duplicate keys in photon heap");
533     r--;
534     break;
535     }
536     }
537    
538     if (heap [heapIdx [r]].pos [d] > p [d])
539     r--;
540     else break;
541     }
542    
543     /* Swap indices (not the nodes they point to) */
544     n2 = heapIdx [l];
545     heapIdx [l] = heapIdx [r];
546     heapIdx [r] = n2;
547     /* Update reverse indices */
548     heapXdi [heapIdx [l]] = l;
549     heapXdi [n2] = r;
550     } while (l < r);
551    
552     /* Swap indices of convergence and pivot nodes */
553     heapIdx [r] = heapIdx [l];
554     heapIdx [l] = heapIdx [right];
555     heapIdx [right] = n2;
556     /* Update reverse indices */
557     heapXdi [heapIdx [r]] = r;
558     heapXdi [heapIdx [l]] = l;
559     heapXdi [n2] = right;
560     if (l >= m) right = l - 1;
561     if (l <= m) left = l + 1;
562     }
563    
564     /* Once left & right have converged at m, we have found the median */
565     return m;
566     }
567    
568    
569    
570     void buildHeap (Photon* heap, unsigned long* heapIdx,
571     unsigned long* heapXdi, const float min [3],
572     const float max [3], unsigned long left,
573     unsigned long right, unsigned long root)
574     /* Recursive part of balancePhotons(..). Builds heap from subarray
575     defined by indices left and right. min and max are the minimum resp.
576     maximum photon positions in the array. root is the index of the
577     current subtree's root, which corresponds to the median's 1-based
578     index in the heap. heapIdx are the balanced heap indices. The heap
579     is accessed indirectly through these. heapXdi are the reverse indices
580     from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */
581     {
582     float maxLeft [3], minRight [3];
583     Photon rootNode;
584     unsigned d;
585    
586     /* Choose median for dimension with largest spread and partition
587     accordingly */
588     const float d0 = max [0] - min [0],
589     d1 = max [1] - min [1],
590     d2 = max [2] - min [2];
591     const unsigned char dim = d0 > d1 ? d0 > d2 ? 0 : 2
592     : d1 > d2 ? 1 : 2;
593     const unsigned long median =
594     left == right ? left
595     : medianPartition(heap, heapIdx, heapXdi, left, right, dim);
596    
597     /* Place median at root of current subtree. This consists of swapping
598     the median and the root nodes and updating the heap indices */
599     memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon));
600     memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon));
601     setPhotonDiscr(rootNode, dim);
602     memcpy(heap + root - 1, &rootNode, sizeof(Photon));
603     heapIdx [heapXdi [root - 1]] = heapIdx [median];
604     heapXdi [heapIdx [median]] = heapXdi [root - 1];
605     heapIdx [median] = root - 1;
606     heapXdi [root - 1] = median;
607    
608     /* Update bounds for left and right subtrees and recurse on them */
609     for (d = 0; d <= 2; d++)
610     if (d == dim)
611     maxLeft [d] = minRight [d] = rootNode.pos [d];
612     else {
613     maxLeft [d] = max [d];
614     minRight [d] = min [d];
615     }
616    
617     if (left < median)
618     buildHeap(heap, heapIdx, heapXdi, min, maxLeft,
619     left, median - 1, root << 1);
620    
621     if (right > median)
622     buildHeap(heap, heapIdx, heapXdi, minRight, max,
623     median + 1, right, (root << 1) + 1);
624     }
625    
626    
627    
628     void balancePhotons (PhotonMap* pmap, double *photonFlux)
629     {
630     Photon *heap = pmap -> heap;
631     unsigned long i;
632     unsigned long *heapIdx; /* Photon index array */
633     unsigned long *heapXdi; /* Reverse index to heapIdx */
634     unsigned j;
635     COLOR flux;
636     /* Need doubles here to reduce errors from increment */
637     double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
638     FVECT d;
639    
640     if (pmap -> heapEnd) {
641     pmap -> heapSize = pmap -> heapEnd;
642     heapIdx = (unsigned long*)malloc(pmap -> heapSize *
643     sizeof(unsigned long));
644     heapXdi = (unsigned long*)malloc(pmap -> heapSize *
645     sizeof(unsigned long));
646     if (!heapIdx || !heapXdi)
647     error(USER, "can't allocate heap index");
648    
649     for (i = 0; i < pmap -> heapSize; i++) {
650     /* Initialize index arrays */
651     heapXdi [i] = heapIdx [i] = i;
652     getPhotonFlux(heap + i, flux);
653    
654     /* Scale photon's flux (hitherto normalised to 1 over RGB); in case
655     * of a contrib photon map, this is done per light source, and
656     * photonFlux is assumed to be an array */
657     if (photonFlux) {
658     scalecolor(flux, photonFlux [isContribPmap(pmap) ?
659     photonSrcIdx(pmap, heap + i) : 0]);
660     setPhotonFlux(heap + i, flux);
661     }
662    
663     /* Need a double here */
664     addcolor(avgFlux, flux);
665    
666     /* Add photon position to centre of gravity */
667     for (j = 0; j < 3; j++)
668     CoG [j] += heap [i].pos [j];
669     }
670    
671     /* Average photon positions to get centre of gravity */
672     for (j = 0; j < 3; j++)
673     pmap -> CoG [j] = CoG [j] /= pmap -> heapSize;
674    
675     /* Compute average photon distance to CoG */
676     for (i = 0; i < pmap -> heapSize; i++) {
677     VSUB(d, heap [i].pos, CoG);
678     CoGdist += DOT(d, d);
679     }
680    
681     pmap -> CoGdist = CoGdist /= pmap -> heapSize;
682    
683     /* Average photon flux based on RGBE representation */
684     scalecolor(avgFlux, 1.0 / pmap -> heapSize);
685     copycolor(pmap -> photonFlux, avgFlux);
686    
687     /* Build kd-tree */
688     buildHeap(pmap -> heap, heapIdx, heapXdi, pmap -> minPos,
689     pmap -> maxPos, 0, pmap -> heapSize - 1, 1);
690    
691     free(heapIdx);
692     free(heapXdi);
693     }
694     }
695    
696    
697    
698     void deletePhotons (PhotonMap* pmap)
699     {
700     free(pmap -> heap);
701     free(pmap -> squeue);
702     free(pmap -> biasCompHist);
703    
704     pmap -> heapSize = 0;
705     pmap -> minGather = pmap -> maxGather =
706     pmap -> squeueSize = pmap -> squeueEnd = 0;
707     }