ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapdata.c
Revision: 2.1
Committed: Tue Feb 24 19:39:26 2015 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial check-in of photon map addition by Roland Schregle

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     $Id: pmapdata.c,v 4.30 2015/01/30 19:22:10 taschreg Exp taschreg $
11     */
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     else {
388     sprintf(errmsg, "max photon lookup radius clamped to %.1e",
389     pmap -> maxDist0);
390     error(WARNING, errmsg);
391     }
392     }
393    
394     /* Reset successful lookup counter */
395     pmap -> numLookups = 0;
396     }
397     else {
398     /* Increment successful lookup counter and reduce max search radius if
399     * wraparound */
400     pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT;
401     if (!pmap -> numLookups)
402     pmap -> maxDist0 *= PMAP_MAXDIST_DEC;
403    
404     redo = 0;
405     }
406     } while (redo);
407     }
408    
409    
410    
411     static void nearest1Neighbour (PhotonMap *pmap, const float pos [3],
412     const float norm [3], Photon **photon,
413     unsigned long node)
414     /* Recursive part of find1Photon(..).
415     Note that all heap index handling is 1-based, but accesses to the
416     arrays are 0-based! */
417     {
418     Photon *p = pmap -> heap + node - 1;
419     /* Signed distance to current photon's splitting plane */
420     float d = pos [photonDiscr(*p)] - p -> pos [photonDiscr(*p)],
421     d2 = d * d;
422     float dv [3];
423    
424     /* Search subtree closer to pos first; exclude other subtree if the
425     distance to the splitting plane is greater than maxDist */
426     if (d < 0) {
427     if (node << 1 <= pmap -> heapSize)
428     nearest1Neighbour(pmap, pos, norm, photon, node << 1);
429     if (d2 < pmap -> maxDist && node << 1 < pmap -> heapSize)
430     nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
431     }
432     else {
433     if (node << 1 < pmap -> heapSize)
434     nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
435     if (d2 < pmap -> maxDist && node << 1 <= pmap -> heapSize)
436     nearest1Neighbour(pmap, pos, norm, photon, node << 1);
437     }
438    
439     /* Squared distance to current photon */
440     dv [0] = pos [0] - p -> pos [0];
441     dv [1] = pos [1] - p -> pos [1];
442     dv [2] = pos [2] - p -> pos [2];
443     d2 = DOT(dv, dv);
444    
445     if (d2 < pmap -> maxDist && DOT(norm, p -> norm) > 0) {
446     /* Closest photon so far with similar normal */
447     pmap -> maxDist = d2;
448     *photon = p;
449     }
450     }
451    
452    
453    
454     Photon* find1Photon (PhotonMap *pmap, const RAY* ray)
455     {
456     float fpos [3], norm [3];
457     Photon* photon = NULL;
458    
459     VCOPY(fpos, ray -> rop);
460     VCOPY(norm, ray -> ron);
461     pmap -> maxDist = thescene.cusize;
462     nearest1Neighbour(pmap, fpos, norm, &photon, 1);
463    
464     return photon;
465     }
466    
467    
468    
469     static unsigned long medianPartition (const Photon* heap,
470     unsigned long* heapIdx,
471     unsigned long* heapXdi,
472     unsigned long left,
473     unsigned long right, unsigned dim)
474     /* Returns index to median in heap from indices left to right
475     (inclusive) in dimension dim. The heap is partitioned relative to
476     median using a quicksort algorithm. The heap indices in heapIdx are
477     sorted rather than the heap itself. */
478     {
479     register const float* p;
480     const unsigned long n = right - left + 1;
481     register unsigned long l, r, lg2, n2, m;
482     register unsigned d;
483    
484     /* Round down n to nearest power of 2 */
485     for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2);
486     n2 = 1 << lg2;
487    
488     /* Determine median position; this takes into account the fact that
489     only the last level in the heap can be partially empty, and that
490     it fills from left to right */
491     m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1));
492    
493     while (right > left) {
494     /* Pivot node */
495     p = heap [heapIdx [right]].pos;
496     l = left;
497     r = right - 1;
498    
499     /* l & r converge, swapping elements out of order with respect to
500     pivot node. Identical keys are resolved by cycling through
501     dim. The convergence point is then the pivot's position. */
502     do {
503     while (l <= r) {
504     d = dim;
505    
506     while (heap [heapIdx [l]].pos [d] == p [d]) {
507     d = (d + 1) % 3;
508    
509     if (d == dim) {
510     /* Ignore dupes? */
511     error(WARNING, "duplicate keys in photon heap");
512     l++;
513     break;
514     }
515     }
516    
517     if (heap [heapIdx [l]].pos [d] < p [d])
518     l++;
519     else break;
520     }
521    
522     while (r > l) {
523     d = dim;
524    
525     while (heap [heapIdx [r]].pos [d] == p [d]) {
526     d = (d + 1) % 3;
527    
528     if (d == dim) {
529     /* Ignore dupes? */
530     error(WARNING, "duplicate keys in photon heap");
531     r--;
532     break;
533     }
534     }
535    
536     if (heap [heapIdx [r]].pos [d] > p [d])
537     r--;
538     else break;
539     }
540    
541     /* Swap indices (not the nodes they point to) */
542     n2 = heapIdx [l];
543     heapIdx [l] = heapIdx [r];
544     heapIdx [r] = n2;
545     /* Update reverse indices */
546     heapXdi [heapIdx [l]] = l;
547     heapXdi [n2] = r;
548     } while (l < r);
549    
550     /* Swap indices of convergence and pivot nodes */
551     heapIdx [r] = heapIdx [l];
552     heapIdx [l] = heapIdx [right];
553     heapIdx [right] = n2;
554     /* Update reverse indices */
555     heapXdi [heapIdx [r]] = r;
556     heapXdi [heapIdx [l]] = l;
557     heapXdi [n2] = right;
558     if (l >= m) right = l - 1;
559     if (l <= m) left = l + 1;
560     }
561    
562     /* Once left & right have converged at m, we have found the median */
563     return m;
564     }
565    
566    
567    
568     void buildHeap (Photon* heap, unsigned long* heapIdx,
569     unsigned long* heapXdi, const float min [3],
570     const float max [3], unsigned long left,
571     unsigned long right, unsigned long root)
572     /* Recursive part of balancePhotons(..). Builds heap from subarray
573     defined by indices left and right. min and max are the minimum resp.
574     maximum photon positions in the array. root is the index of the
575     current subtree's root, which corresponds to the median's 1-based
576     index in the heap. heapIdx are the balanced heap indices. The heap
577     is accessed indirectly through these. heapXdi are the reverse indices
578     from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */
579     {
580     float maxLeft [3], minRight [3];
581     Photon rootNode;
582     unsigned d;
583    
584     /* Choose median for dimension with largest spread and partition
585     accordingly */
586     const float d0 = max [0] - min [0],
587     d1 = max [1] - min [1],
588     d2 = max [2] - min [2];
589     const unsigned char dim = d0 > d1 ? d0 > d2 ? 0 : 2
590     : d1 > d2 ? 1 : 2;
591     const unsigned long median =
592     left == right ? left
593     : medianPartition(heap, heapIdx, heapXdi, left, right, dim);
594    
595     /* Place median at root of current subtree. This consists of swapping
596     the median and the root nodes and updating the heap indices */
597     memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon));
598     memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon));
599     setPhotonDiscr(rootNode, dim);
600     memcpy(heap + root - 1, &rootNode, sizeof(Photon));
601     heapIdx [heapXdi [root - 1]] = heapIdx [median];
602     heapXdi [heapIdx [median]] = heapXdi [root - 1];
603     heapIdx [median] = root - 1;
604     heapXdi [root - 1] = median;
605    
606     /* Update bounds for left and right subtrees and recurse on them */
607     for (d = 0; d <= 2; d++)
608     if (d == dim)
609     maxLeft [d] = minRight [d] = rootNode.pos [d];
610     else {
611     maxLeft [d] = max [d];
612     minRight [d] = min [d];
613     }
614    
615     if (left < median)
616     buildHeap(heap, heapIdx, heapXdi, min, maxLeft,
617     left, median - 1, root << 1);
618    
619     if (right > median)
620     buildHeap(heap, heapIdx, heapXdi, minRight, max,
621     median + 1, right, (root << 1) + 1);
622     }
623    
624    
625    
626     void balancePhotons (PhotonMap* pmap, double *photonFlux)
627     {
628     Photon *heap = pmap -> heap;
629     unsigned long i;
630     unsigned long *heapIdx; /* Photon index array */
631     unsigned long *heapXdi; /* Reverse index to heapIdx */
632     unsigned j;
633     COLOR flux;
634     /* Need doubles here to reduce errors from increment */
635     double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
636     FVECT d;
637    
638     if (pmap -> heapEnd) {
639     pmap -> heapSize = pmap -> heapEnd;
640     heapIdx = (unsigned long*)malloc(pmap -> heapSize *
641     sizeof(unsigned long));
642     heapXdi = (unsigned long*)malloc(pmap -> heapSize *
643     sizeof(unsigned long));
644     if (!heapIdx || !heapXdi)
645     error(USER, "can't allocate heap index");
646    
647     for (i = 0; i < pmap -> heapSize; i++) {
648     /* Initialize index arrays */
649     heapXdi [i] = heapIdx [i] = i;
650     getPhotonFlux(heap + i, flux);
651    
652     /* Scale photon's flux (hitherto normalised to 1 over RGB); in case
653     * of a contrib photon map, this is done per light source, and
654     * photonFlux is assumed to be an array */
655     if (photonFlux) {
656     scalecolor(flux, photonFlux [isContribPmap(pmap) ?
657     photonSrcIdx(pmap, heap + i) : 0]);
658     setPhotonFlux(heap + i, flux);
659     }
660    
661     /* Need a double here */
662     addcolor(avgFlux, flux);
663    
664     /* Add photon position to centre of gravity */
665     for (j = 0; j < 3; j++)
666     CoG [j] += heap [i].pos [j];
667     }
668    
669     /* Average photon positions to get centre of gravity */
670     for (j = 0; j < 3; j++)
671     pmap -> CoG [j] = CoG [j] /= pmap -> heapSize;
672    
673     /* Compute average photon distance to CoG */
674     for (i = 0; i < pmap -> heapSize; i++) {
675     VSUB(d, heap [i].pos, CoG);
676     CoGdist += DOT(d, d);
677     }
678    
679     pmap -> CoGdist = CoGdist /= pmap -> heapSize;
680    
681     /* Average photon flux based on RGBE representation */
682     scalecolor(avgFlux, 1.0 / pmap -> heapSize);
683     copycolor(pmap -> photonFlux, avgFlux);
684    
685     /* Build kd-tree */
686     buildHeap(pmap -> heap, heapIdx, heapXdi, pmap -> minPos,
687     pmap -> maxPos, 0, pmap -> heapSize - 1, 1);
688    
689     free(heapIdx);
690     free(heapXdi);
691     }
692     }
693    
694    
695    
696     void deletePhotons (PhotonMap* pmap)
697     {
698     free(pmap -> heap);
699     free(pmap -> squeue);
700     free(pmap -> biasCompHist);
701    
702     pmap -> heapSize = 0;
703     pmap -> minGather = pmap -> maxGather =
704     pmap -> squeueSize = pmap -> squeueEnd = 0;
705     }