ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapkdt.c
Revision: 1.2
Committed: Mon Aug 14 21:12:10 2017 UTC (6 years, 8 months ago) by rschregle
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R1
Changes since 1.1: +7 -9 lines
Log Message:
Updated photon map code for Windows; no multproc or ooC for now

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.1 2016/05/18 08:22:45 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) && pmap -> srcContrib) {
350 /* Lookup in contribution photon map */
351 OBJREC *srcMod;
352 const int srcIdx = photonSrcIdx(pmap, p);
353
354 if (srcIdx < 0 || srcIdx >= nsources)
355 error(INTERNAL, "invalid light source index in photon map");
356
357 srcMod = findmaterial(source [srcIdx].so);
358
359 /* Reject photon if contributions from light source which emitted it
360 * are not sought */
361 if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
362 return;
363
364 /* Reject non-caustic photon if lookup for caustic contribs */
365 if (pmap -> lookupCaustic & !p -> caustic)
366 return;
367 }
368
369 /* Squared distance to current photon (note dist2() requires doubles) */
370 VSUB(dv, pos, p -> pos);
371 d2 = DOT(dv, dv);
372
373 /* Accept photon if closer than current max dist & add to priority queue */
374 if (d2 < pmap -> maxDist2) {
375 if (pmap -> squeue.tail < sqSize) {
376 /* Priority queue not full; append photon and restore heap */
377 i = ++pmap -> squeue.tail;
378
379 while (i > 1 && sq [(i >> 1) - 1].dist2 <= d2) {
380 sq [i - 1].idx = sq [(i >> 1) - 1].idx;
381 sq [i - 1].dist2 = sq [(i >> 1) - 1].dist2;
382 i >>= 1;
383 }
384
385 sq [--i].idx = (PhotonIdx)p;
386 sq [i].dist2 = d2;
387 /* Update maxDist if we've just filled the queue */
388 if (pmap -> squeue.tail >= pmap -> squeue.len)
389 pmap -> maxDist2 = sq [0].dist2;
390 }
391 else {
392 /* Priority queue full; replace maximum, restore heap, and
393 update maxDist */
394 i = 1;
395
396 while (i <= sqSize >> 1) {
397 j = i << 1;
398 if (j < sqSize && sq [j - 1].dist2 < sq [j].dist2)
399 j++;
400 if (d2 >= sq [j - 1].dist2)
401 break;
402 sq [i - 1].idx = sq [j - 1].idx;
403 sq [i - 1].dist2 = sq [j - 1].dist2;
404 i = j;
405 }
406
407 sq [--i].idx = (PhotonIdx)p;
408 sq [i].dist2 = d2;
409 pmap -> maxDist2 = sq [0].dist2;
410 }
411 }
412 }
413
414
415
416 void kdT_FindPhotons (struct PhotonMap *pmap, const FVECT pos,
417 const FVECT norm)
418 {
419 float p [3], n [3];
420
421 /* Photon pos & normal stored at lower precision */
422 VCOPY(p, pos);
423 VCOPY(n, norm);
424 kdT_FindNearest(pmap, p, n, 1);
425 }
426
427
428
429 static void kdT_Find1Nearest (PhotonMap *pmap, const float pos [3],
430 const float norm [3], Photon **photon,
431 unsigned long node)
432 /* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to
433 * pos with similar normal. Note that all heap and queue indices are
434 * 1-based, but accesses to the arrays are 0-based! */
435 {
436 Photon *p = (Photon*)pmap -> store.nodes + node - 1;
437 /* Signed distance to current photon's splitting plane */
438 float d = pos [p -> discr] - p -> pos [p -> discr], d2 = d * d,
439 dv [3];
440
441 /* Search subtree closer to pos first; exclude other subtree if the
442 distance to the splitting plane is greater than maxDist */
443 if (d < 0) {
444 if (node << 1 <= pmap -> numPhotons)
445 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
446
447 if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)
448 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
449 }
450 else {
451 if (node << 1 < pmap -> numPhotons)
452 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
453
454 if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)
455 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
456 }
457
458 /* Squared distance to current photon */
459 VSUB(dv, pos, p -> pos);
460 d2 = DOT(dv, dv);
461
462 if (d2 < pmap -> maxDist2 &&
463 DOT(norm, p -> norm) > PMAP_NORM_TOL * 127 * frandom()) {
464 /* Closest photon so far with similar normal. We allow for tolerance
465 * to account for perturbation in the latter; note the photon normal
466 * is coded in the range [-127,127], hence we factor this in */
467 pmap -> maxDist2 = d2;
468 *photon = p;
469 }
470 }
471
472
473
474 void kdT_Find1Photon (struct PhotonMap *pmap, const FVECT pos,
475 const FVECT norm, Photon *photon)
476 {
477 float p [3], n [3];
478 Photon *pnn;
479
480 /* Photon pos & normal stored at lower precision */
481 VCOPY(p, pos);
482 VCOPY(n, norm);
483 kdT_Find1Nearest(pmap, p, n, &pnn, 1);
484 memcpy(photon, pnn, sizeof(Photon));
485 }
486
487
488
489 int kdT_GetPhoton (const struct PhotonMap *pmap, PhotonIdx idx,
490 Photon *photon)
491 {
492 memcpy(photon, idx, sizeof(Photon));
493 return 0;
494 }
495
496
497
498 Photon *kdT_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
499 {
500 return idx;
501 }
502
503
504
505 PhotonIdx kdT_FirstPhoton (const struct PhotonMap* pmap)
506 {
507 return pmap -> store.nodes;
508 }
509
510
511
512 void kdT_Delete (PhotonKdTree *kdt)
513 {
514 free(kdt -> nodes);
515 kdt -> nodes = NULL;
516 }