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 (9 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial check-in of photon map addition by Roland Schregle

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