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

# Content
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.31 2015/04/22 15:37:07 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 #ifdef PMAP_LOOKUP_REDO
388 else {
389 sprintf(errmsg, "max photon lookup radius clamped to %.1e",
390 pmap -> maxDist0);
391 error(WARNING, errmsg);
392 }
393 #endif
394 }
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 }