ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapkdt.c
Revision: 1.3
Committed: Wed Jan 24 19:39:05 2018 UTC (6 years, 3 months ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 1.2: +19 -15 lines
Log Message:
Hack to render contrib pmap as regular global pmap with rpict/rtrace/rvu

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.2 2017/08/14 21:12:10 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 VCOPY(n, norm);
428 kdT_FindNearest(pmap, p, n, 1);
429 }
430
431
432
433 static void kdT_Find1Nearest (PhotonMap *pmap, const float pos [3],
434 const float norm [3], Photon **photon,
435 unsigned long node)
436 /* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to
437 * pos with similar normal. Note that all heap and queue indices are
438 * 1-based, but accesses to the arrays are 0-based! */
439 {
440 Photon *p = (Photon*)pmap -> store.nodes + node - 1;
441 /* Signed distance to current photon's splitting plane */
442 float d = pos [p -> discr] - p -> pos [p -> discr], d2 = d * d,
443 dv [3];
444
445 /* Search subtree closer to pos first; exclude other subtree if the
446 distance to the splitting plane is greater than maxDist */
447 if (d < 0) {
448 if (node << 1 <= pmap -> numPhotons)
449 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
450
451 if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)
452 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
453 }
454 else {
455 if (node << 1 < pmap -> numPhotons)
456 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
457
458 if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)
459 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
460 }
461
462 /* Squared distance to current photon */
463 VSUB(dv, pos, p -> pos);
464 d2 = DOT(dv, dv);
465
466 if (d2 < pmap -> maxDist2 &&
467 DOT(norm, p -> norm) > PMAP_NORM_TOL * 127 * frandom()) {
468 /* Closest photon so far with similar normal. We allow for tolerance
469 * to account for perturbation in the latter; note the photon normal
470 * is coded in the range [-127,127], hence we factor this in */
471 pmap -> maxDist2 = d2;
472 *photon = p;
473 }
474 }
475
476
477
478 void kdT_Find1Photon (struct PhotonMap *pmap, const FVECT pos,
479 const FVECT norm, Photon *photon)
480 {
481 float p [3], n [3];
482 Photon *pnn;
483
484 /* Photon pos & normal stored at lower precision */
485 VCOPY(p, pos);
486 VCOPY(n, norm);
487 kdT_Find1Nearest(pmap, p, n, &pnn, 1);
488 memcpy(photon, pnn, sizeof(Photon));
489 }
490
491
492
493 int kdT_GetPhoton (const struct PhotonMap *pmap, PhotonIdx idx,
494 Photon *photon)
495 {
496 memcpy(photon, idx, sizeof(Photon));
497 return 0;
498 }
499
500
501
502 Photon *kdT_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
503 {
504 return idx;
505 }
506
507
508
509 PhotonIdx kdT_FirstPhoton (const struct PhotonMap* pmap)
510 {
511 return pmap -> store.nodes;
512 }
513
514
515
516 void kdT_Delete (PhotonKdTree *kdt)
517 {
518 free(kdt -> nodes);
519 kdt -> nodes = NULL;
520 }