ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapkdt.c
Revision: 1.6
Committed: Wed Apr 8 15:14:21 2020 UTC (4 years ago) by rschregle
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R3, HEAD
Changes since 1.5: +18 -8 lines
Log Message:
Fixed est00pid bug in single photon lookups that returned junk when none found, added code to detect and handle (ignore).

File Contents

# Content
1 /*
2 ======================================================================
3 In-core kd-tree for photon map
4
5 Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
6 (c) Fraunhofer Institute for Solar Energy Systems,
7 (c) Lucerne University of Applied Sciences and Arts,
8 supported by the Swiss National Science Foundation (SNSF, #147053)
9 ======================================================================
10
11 $Id: pmapkdt.c,v 1.5 2018/11/08 00:54:07 greg Exp $
12 */
13
14
15
16 #include "pmapdata.h" /* Includes pmapkdt.h */
17 #include "source.h"
18 #include "otspecial.h"
19 #include "random.h"
20
21
22
23
24 void kdT_Null (PhotonKdTree *kdt)
25 {
26 kdt -> nodes = NULL;
27 }
28
29
30
31 static unsigned long kdT_MedianPartition (const Photon *heap,
32 unsigned long *heapIdx,
33 unsigned long *heapXdi,
34 unsigned long left,
35 unsigned long right, unsigned dim)
36 /* Returns index to median in heap from indices left to right
37 (inclusive) in dimension dim. The heap is partitioned relative to
38 median using a quicksort algorithm. The heap indices in heapIdx are
39 sorted rather than the heap itself. */
40 {
41 const float *p;
42 unsigned long l, r, lg2, n2, m, n = right - left + 1;
43 unsigned d;
44
45 /* Round down n to nearest power of 2 */
46 for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2);
47 n2 = 1 << lg2;
48
49 /* Determine median position; this takes into account the fact that
50 only the last level in the heap can be partially empty, and that
51 it fills from left to right */
52 m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1));
53
54 while (right > left) {
55 /* Pivot node */
56 p = heap [heapIdx [right]].pos;
57 l = left;
58 r = right - 1;
59
60 /* l & r converge, swapping elements out of order with respect to
61 pivot node. Identical keys are resolved by cycling through
62 dim. The convergence point is then the pivot's position. */
63 do {
64 while (l <= r) {
65 d = dim;
66
67 while (heap [heapIdx [l]].pos [d] == p [d]) {
68 d = (d + 1) % 3;
69
70 if (d == dim) {
71 /* Ignore dupes? */
72 error(WARNING, "duplicate keys in photon heap");
73 l++;
74 break;
75 }
76 }
77
78 if (heap [heapIdx [l]].pos [d] < p [d])
79 l++;
80 else break;
81 }
82
83 while (r > l) {
84 d = dim;
85
86 while (heap [heapIdx [r]].pos [d] == p [d]) {
87 d = (d + 1) % 3;
88
89 if (d == dim) {
90 /* Ignore dupes? */
91 error(WARNING, "duplicate keys in photon heap");
92 r--;
93 break;
94 }
95 }
96
97 if (heap [heapIdx [r]].pos [d] > p [d])
98 r--;
99 else break;
100 }
101
102 /* Swap indices (not the nodes they point to) */
103 n2 = heapIdx [l];
104 heapIdx [l] = heapIdx [r];
105 heapIdx [r] = n2;
106 /* Update reverse indices */
107 heapXdi [heapIdx [l]] = l;
108 heapXdi [n2] = r;
109 } while (l < r);
110
111 /* Swap indices of convergence and pivot nodes */
112 heapIdx [r] = heapIdx [l];
113 heapIdx [l] = heapIdx [right];
114 heapIdx [right] = n2;
115 /* Update reverse indices */
116 heapXdi [heapIdx [r]] = r;
117 heapXdi [heapIdx [l]] = l;
118 heapXdi [n2] = right;
119
120 if (l >= m)
121 right = l - 1;
122 if (l <= m)
123 left = l + 1;
124 }
125
126 /* Once left & right have converged at m, we have found the median */
127 return m;
128 }
129
130
131
132 static void kdT_Build (Photon *heap, unsigned long *heapIdx,
133 unsigned long *heapXdi, const float min [3],
134 const float max [3], unsigned long left,
135 unsigned long right, unsigned long root)
136 /* Recursive part of balancePhotons(..). Builds heap from subarray
137 defined by indices left and right. min and max are the minimum resp.
138 maximum photon positions in the array. root is the index of the
139 current subtree's root, which corresponds to the median's 1-based
140 index in the heap. heapIdx are the balanced heap indices. The heap
141 is accessed indirectly through these. heapXdi are the reverse indices
142 from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */
143 {
144 float maxLeft [3], minRight [3];
145 Photon rootNode;
146 unsigned d;
147
148 /* Choose median for dimension with largest spread and partition
149 accordingly */
150 const float d0 = max [0] - min [0],
151 d1 = max [1] - min [1],
152 d2 = max [2] - min [2];
153 const unsigned char dim = d0 > d1 ? d0 > d2 ? 0 : 2
154 : d1 > d2 ? 1 : 2;
155 const unsigned long median = left == right
156 ? left
157 : kdT_MedianPartition(heap, heapIdx, heapXdi,
158 left, right, dim);
159
160 /* Place median at root of current subtree. This consists of swapping
161 the median and the root nodes and updating the heap indices */
162 memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon));
163 memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon));
164 rootNode.discr = dim;
165 memcpy(heap + root - 1, &rootNode, sizeof(Photon));
166 heapIdx [heapXdi [root - 1]] = heapIdx [median];
167 heapXdi [heapIdx [median]] = heapXdi [root - 1];
168 heapIdx [median] = root - 1;
169 heapXdi [root - 1] = median;
170
171 /* Update bounds for left and right subtrees and recurse on them */
172 for (d = 0; d <= 2; d++)
173 if (d == dim)
174 maxLeft [d] = minRight [d] = rootNode.pos [d];
175 else {
176 maxLeft [d] = max [d];
177 minRight [d] = min [d];
178 }
179
180 if (left < median)
181 kdT_Build(heap, heapIdx, heapXdi, min, maxLeft, left, median - 1,
182 root << 1);
183
184 if (right > median)
185 kdT_Build(heap, heapIdx, heapXdi, minRight, max, median + 1, right,
186 (root << 1) + 1);
187 }
188
189
190
191 void kdT_BuildPhotonMap (struct PhotonMap *pmap)
192 {
193 Photon *nodes;
194 unsigned long i;
195 unsigned long *heapIdx, /* Photon index array */
196 *heapXdi; /* Reverse index to heapIdx */
197
198 /* Allocate kd-tree nodes and load photons from heap file */
199 if (!(nodes = calloc(pmap -> numPhotons, sizeof(Photon))))
200 error(SYSTEM, "failed in-core heap allocation in kdT_BuildPhotonMap");
201
202 rewind(pmap -> heap);
203 i = fread(nodes, sizeof(Photon), pmap -> numPhotons, pmap -> heap);
204 if (i !=
205 pmap -> numPhotons)
206 error(SYSTEM, "failed loading photon heap in kdT_BuildPhotonMap");
207
208 pmap -> store.nodes = nodes;
209 heapIdx = calloc(pmap -> numPhotons, sizeof(unsigned long));
210 heapXdi = calloc(pmap -> numPhotons, sizeof(unsigned long));
211 if (!heapIdx || !heapXdi)
212 error(SYSTEM, "failed heap index allocation in kdT_BuildPhotonMap");
213
214 /* Initialize index arrays */
215 for (i = 0; i < pmap -> numPhotons; i++)
216 heapXdi [i] = heapIdx [i] = i;
217
218 /* Build kd-tree */
219 kdT_Build(nodes, heapIdx, heapXdi, pmap -> minPos, pmap -> maxPos, 0,
220 pmap -> numPhotons - 1, 1);
221
222 /* Cleanup */
223 free(heapIdx);
224 free(heapXdi);
225 }
226
227
228
229 int kdT_SavePhotons (const struct PhotonMap *pmap, FILE *out)
230 {
231 unsigned long i, j;
232 Photon *p = (Photon*)pmap -> store.nodes;
233
234 for (i = 0; i < pmap -> numPhotons; i++, p++) {
235 /* Write photon attributes */
236 for (j = 0; j < 3; j++)
237 putflt(p -> pos [j], out);
238
239 /* Bytewise dump otherwise we have portability probs */
240 for (j = 0; j < 3; j++)
241 putint(p -> norm [j], 1, out);
242
243 #ifdef PMAP_FLOAT_FLUX
244 for (j = 0; j < 3; j++)
245 putflt(p -> flux [j], out);
246 #else
247 for (j = 0; j < 4; j++)
248 putint(p -> flux [j], 1, out);
249 #endif
250
251 putint(p -> primary, sizeof(p -> primary), out);
252 putint(p -> flags, 1, out);
253
254 if (ferror(out))
255 return -1;
256 }
257
258 return 0;
259 }
260
261
262
263 int kdT_LoadPhotons (struct PhotonMap *pmap, FILE *in)
264 {
265 unsigned long i, j;
266 Photon *p;
267
268 /* Allocate kd-tree based on initialised pmap -> numPhotons */
269 pmap -> store.nodes = calloc(sizeof(Photon), pmap -> numPhotons);
270 if (!pmap -> store.nodes)
271 error(SYSTEM, "failed kd-tree allocation in kdT_LoadPhotons");
272
273 /* Get photon attributes */
274 for (i = 0, p = pmap -> store.nodes; i < pmap -> numPhotons; i++, p++) {
275 for (j = 0; j < 3; j++)
276 p -> pos [j] = getflt(in);
277
278 /* Bytewise grab otherwise we have portability probs */
279 for (j = 0; j < 3; j++)
280 p -> norm [j] = getint(1, in);
281
282 #ifdef PMAP_FLOAT_FLUX
283 for (j = 0; j < 3; j++)
284 p -> flux [j] = getflt(in);
285 #else
286 for (j = 0; j < 4; j++)
287 p -> flux [j] = getint(1, in);
288 #endif
289
290 p -> primary = getint(sizeof(p -> primary), in);
291 p -> flags = getint(1, in);
292
293 if (feof(in))
294 return -1;
295 }
296
297 return 0;
298 }
299
300
301
302 void kdT_InitFindPhotons (struct PhotonMap *pmap)
303 {
304 pmap -> squeue.len = pmap -> maxGather + 1;
305 pmap -> squeue.node = calloc(pmap -> squeue.len,
306 sizeof(PhotonSearchQueueNode));
307 if (!pmap -> squeue.node)
308 error(SYSTEM, "can't allocate photon search queue");
309 }
310
311
312
313 static void kdT_FindNearest (PhotonMap *pmap, const float pos [3],
314 const float norm [3], unsigned long node)
315 /* Recursive part of kdT_FindPhotons(). Locate pmap -> squeue.len nearest
316 * neighbours to pos with similar normal and return in search queue starting
317 * at pmap -> squeue.node. Note that all heap and queue indices are
318 * 1-based, but accesses to the arrays are 0-based! */
319 {
320 Photon *p = (Photon*)pmap -> store.nodes + node - 1;
321 unsigned i, j;
322 /* Signed distance to current photon's splitting plane */
323 float d = pos [p -> discr] - p -> pos [p -> discr],
324 d2 = d * d, dv [3];
325 PhotonSearchQueueNode* sq = pmap -> squeue.node;
326 const unsigned sqSize = pmap -> squeue.len;
327
328 /* Search subtree closer to pos first; exclude other subtree if the
329 distance to the splitting plane is greater than maxDist */
330 if (d < 0) {
331 if (node << 1 <= pmap -> numPhotons)
332 kdT_FindNearest(pmap, pos, norm, node << 1);
333
334 if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)
335 kdT_FindNearest(pmap, pos, norm, (node << 1) + 1);
336 }
337 else {
338 if (node << 1 < pmap -> numPhotons)
339 kdT_FindNearest(pmap, pos, norm, (node << 1) + 1);
340
341 if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)
342 kdT_FindNearest(pmap, pos, norm, node << 1);
343 }
344
345 /* Reject photon if normal faces away (ignored for volume photons) with
346 * tolerance to account for perturbation; note photon normal is coded
347 * in range [-127,127], hence we factor this in */
348 if (norm && DOT(norm, p -> norm) <= PMAP_NORM_TOL * 127 * frandom())
349 return;
350
351 if (isContribPmap(pmap)) {
352 /* Lookup in contribution photon map; filter according to emitting
353 * light source if contrib list set, else accept all */
354
355 if (pmap -> srcContrib) {
356 OBJREC *srcMod;
357 const int srcIdx = photonSrcIdx(pmap, p);
358
359 if (srcIdx < 0 || srcIdx >= nsources)
360 error(INTERNAL, "invalid light source index in photon map");
361
362 srcMod = findmaterial(source [srcIdx].so);
363
364 /* Reject photon if contributions from light source which emitted it
365 * are not sought */
366 if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
367 return;
368 }
369
370 /* Reject non-caustic photon if lookup for caustic contribs */
371 if (pmap -> lookupCaustic & !p -> caustic)
372 return;
373 }
374
375 /* Squared distance to current photon (note dist2() requires doubles) */
376 VSUB(dv, pos, p -> pos);
377 d2 = DOT(dv, dv);
378
379 /* Accept photon if closer than current max dist & add to priority queue */
380 if (d2 < pmap -> maxDist2) {
381 if (pmap -> squeue.tail < sqSize) {
382 /* Priority queue not full; append photon and restore heap */
383 i = ++pmap -> squeue.tail;
384
385 while (i > 1 && sq [(i >> 1) - 1].dist2 <= d2) {
386 sq [i - 1].idx = sq [(i >> 1) - 1].idx;
387 sq [i - 1].dist2 = sq [(i >> 1) - 1].dist2;
388 i >>= 1;
389 }
390
391 sq [--i].idx = (PhotonIdx)p;
392 sq [i].dist2 = d2;
393 /* Update maxDist if we've just filled the queue */
394 if (pmap -> squeue.tail >= pmap -> squeue.len)
395 pmap -> maxDist2 = sq [0].dist2;
396 }
397 else {
398 /* Priority queue full; replace maximum, restore heap, and
399 update maxDist */
400 i = 1;
401
402 while (i <= sqSize >> 1) {
403 j = i << 1;
404 if (j < sqSize && sq [j - 1].dist2 < sq [j].dist2)
405 j++;
406 if (d2 >= sq [j - 1].dist2)
407 break;
408 sq [i - 1].idx = sq [j - 1].idx;
409 sq [i - 1].dist2 = sq [j - 1].dist2;
410 i = j;
411 }
412
413 sq [--i].idx = (PhotonIdx)p;
414 sq [i].dist2 = d2;
415 pmap -> maxDist2 = sq [0].dist2;
416 }
417 }
418 }
419
420
421
422 int kdT_FindPhotons (struct PhotonMap *pmap, const FVECT pos,
423 const FVECT norm)
424 {
425 float p [3], n [3];
426
427 /* Photon pos & normal stored at lower precision */
428 VCOPY(p, pos);
429 if (norm)
430 VCOPY(n, norm);
431 kdT_FindNearest(pmap, p, norm ? n : NULL, 1);
432
433 /* Return success or failure (empty queue => none found) */
434 return pmap -> squeue.tail ? 0 : -1;
435 }
436
437
438
439 static void kdT_Find1Nearest (PhotonMap *pmap, const float pos [3],
440 const float norm [3], Photon **photon,
441 unsigned long node)
442 /* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to
443 * pos with similar normal. Note that all heap and queue indices are
444 * 1-based, but accesses to the arrays are 0-based! */
445 {
446 Photon *p = (Photon*)pmap -> store.nodes + node - 1;
447 /* Signed distance to current photon's splitting plane */
448 float d = pos [p -> discr] - p -> pos [p -> discr], d2 = d * d,
449 dv [3];
450
451 /* Search subtree closer to pos first; exclude other subtree if the
452 distance to the splitting plane is greater than maxDist */
453 if (d < 0) {
454 if (node << 1 <= pmap -> numPhotons)
455 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
456
457 if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)
458 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
459 }
460 else {
461 if (node << 1 < pmap -> numPhotons)
462 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
463
464 if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)
465 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
466 }
467
468 /* Squared distance to current photon */
469 VSUB(dv, pos, p -> pos);
470 d2 = DOT(dv, dv);
471
472 if (d2 < pmap -> maxDist2 &&
473 (!norm || DOT(norm, p -> norm) > PMAP_NORM_TOL * 127 * frandom())) {
474 /* Closest photon so far with similar normal. We allow for tolerance
475 * to account for perturbation in the latter; note the photon normal
476 * is coded in the range [-127,127], hence we factor this in */
477 pmap -> maxDist2 = d2;
478 *photon = p;
479 }
480 }
481
482
483
484 int kdT_Find1Photon (struct PhotonMap *pmap, const FVECT pos,
485 const FVECT norm, Photon *photon)
486 {
487 float p [3], n [3];
488 Photon *pnn = NULL;
489
490 /* Photon pos & normal stored at lower precision */
491 VCOPY(p, pos);
492 if (norm)
493 VCOPY(n, norm);
494 kdT_Find1Nearest(pmap, p, norm ? n : NULL, &pnn, 1);
495 if (!pnn)
496 /* No photon found => failed */
497 return -1;
498 else {
499 /* Copy found photon => successs */
500 memcpy(photon, pnn, sizeof(Photon));
501 return 0;
502 }
503 }
504
505
506
507 int kdT_GetPhoton (const struct PhotonMap *pmap, PhotonIdx idx,
508 Photon *photon)
509 {
510 memcpy(photon, idx, sizeof(Photon));
511 return 0;
512 }
513
514
515
516 Photon *kdT_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
517 {
518 return idx;
519 }
520
521
522
523 PhotonIdx kdT_FirstPhoton (const struct PhotonMap* pmap)
524 {
525 return pmap -> store.nodes;
526 }
527
528
529
530 void kdT_Delete (PhotonKdTree *kdt)
531 {
532 free(kdt -> nodes);
533 kdt -> nodes = NULL;
534 }