ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapkdt.c
Revision: 1.4
Committed: Thu May 31 12:34:16 2018 UTC (6 years ago) by rschregle
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2
Changes since 1.3: +8 -6 lines
Log Message:
Fixed volume photon lookup bug (dereferenced NULL normal)

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