ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapkdt.c
Revision: 1.1
Committed: Wed May 18 08:22:45 2016 UTC (7 years, 11 months ago) by rschregle
Content type: text/plain
Branch: MAIN
Log Message:
Added missing sources of iC and ooC pmap

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 2016/02/04 20:23:27 taschreg Exp taschreg $
12 */
13
14
15 #ifdef PMAP_OOC
16 /* Checked in pmapdata.h */
17 #undef PMAP_OOC
18 #endif
19 #include "pmapdata.h" /* Includes pmapkdt.h */
20 #include "source.h"
21
22
23
24
25 void kdT_Null (PhotonKdTree *kdt)
26 {
27 kdt -> nodes = NULL;
28 }
29
30
31
32 static unsigned long kdT_MedianPartition (const Photon *heap,
33 unsigned long *heapIdx,
34 unsigned long *heapXdi,
35 unsigned long left,
36 unsigned long right, unsigned dim)
37 /* Returns index to median in heap from indices left to right
38 (inclusive) in dimension dim. The heap is partitioned relative to
39 median using a quicksort algorithm. The heap indices in heapIdx are
40 sorted rather than the heap itself. */
41 {
42 const float *p;
43 unsigned long l, r, lg2, n2, m, n = right - left + 1;
44 unsigned d;
45
46 /* Round down n to nearest power of 2 */
47 for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2);
48 n2 = 1 << lg2;
49
50 /* Determine median position; this takes into account the fact that
51 only the last level in the heap can be partially empty, and that
52 it fills from left to right */
53 m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1));
54
55 while (right > left) {
56 /* Pivot node */
57 p = heap [heapIdx [right]].pos;
58 l = left;
59 r = right - 1;
60
61 /* l & r converge, swapping elements out of order with respect to
62 pivot node. Identical keys are resolved by cycling through
63 dim. The convergence point is then the pivot's position. */
64 do {
65 while (l <= r) {
66 d = dim;
67
68 while (heap [heapIdx [l]].pos [d] == p [d]) {
69 d = (d + 1) % 3;
70
71 if (d == dim) {
72 /* Ignore dupes? */
73 error(WARNING, "duplicate keys in photon heap");
74 l++;
75 break;
76 }
77 }
78
79 if (heap [heapIdx [l]].pos [d] < p [d])
80 l++;
81 else break;
82 }
83
84 while (r > l) {
85 d = dim;
86
87 while (heap [heapIdx [r]].pos [d] == p [d]) {
88 d = (d + 1) % 3;
89
90 if (d == dim) {
91 /* Ignore dupes? */
92 error(WARNING, "duplicate keys in photon heap");
93 r--;
94 break;
95 }
96 }
97
98 if (heap [heapIdx [r]].pos [d] > p [d])
99 r--;
100 else break;
101 }
102
103 /* Swap indices (not the nodes they point to) */
104 n2 = heapIdx [l];
105 heapIdx [l] = heapIdx [r];
106 heapIdx [r] = n2;
107 /* Update reverse indices */
108 heapXdi [heapIdx [l]] = l;
109 heapXdi [n2] = r;
110 } while (l < r);
111
112 /* Swap indices of convergence and pivot nodes */
113 heapIdx [r] = heapIdx [l];
114 heapIdx [l] = heapIdx [right];
115 heapIdx [right] = n2;
116 /* Update reverse indices */
117 heapXdi [heapIdx [r]] = r;
118 heapXdi [heapIdx [l]] = l;
119 heapXdi [n2] = right;
120
121 if (l >= m)
122 right = l - 1;
123 if (l <= m)
124 left = l + 1;
125 }
126
127 /* Once left & right have converged at m, we have found the median */
128 return m;
129 }
130
131
132
133 static void kdT_Build (Photon *heap, unsigned long *heapIdx,
134 unsigned long *heapXdi, const float min [3],
135 const float max [3], unsigned long left,
136 unsigned long right, unsigned long root)
137 /* Recursive part of balancePhotons(..). Builds heap from subarray
138 defined by indices left and right. min and max are the minimum resp.
139 maximum photon positions in the array. root is the index of the
140 current subtree's root, which corresponds to the median's 1-based
141 index in the heap. heapIdx are the balanced heap indices. The heap
142 is accessed indirectly through these. heapXdi are the reverse indices
143 from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */
144 {
145 float maxLeft [3], minRight [3];
146 Photon rootNode;
147 unsigned d;
148
149 /* Choose median for dimension with largest spread and partition
150 accordingly */
151 const float d0 = max [0] - min [0],
152 d1 = max [1] - min [1],
153 d2 = max [2] - min [2];
154 const unsigned char dim = d0 > d1 ? d0 > d2 ? 0 : 2
155 : d1 > d2 ? 1 : 2;
156 const unsigned long median = left == right
157 ? left
158 : kdT_MedianPartition(heap, heapIdx, heapXdi,
159 left, right, dim);
160
161 /* Place median at root of current subtree. This consists of swapping
162 the median and the root nodes and updating the heap indices */
163 memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon));
164 memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon));
165 rootNode.discr = dim;
166 memcpy(heap + root - 1, &rootNode, sizeof(Photon));
167 heapIdx [heapXdi [root - 1]] = heapIdx [median];
168 heapXdi [heapIdx [median]] = heapXdi [root - 1];
169 heapIdx [median] = root - 1;
170 heapXdi [root - 1] = median;
171
172 /* Update bounds for left and right subtrees and recurse on them */
173 for (d = 0; d <= 2; d++)
174 if (d == dim)
175 maxLeft [d] = minRight [d] = rootNode.pos [d];
176 else {
177 maxLeft [d] = max [d];
178 minRight [d] = min [d];
179 }
180
181 if (left < median)
182 kdT_Build(heap, heapIdx, heapXdi, min, maxLeft, left, median - 1,
183 root << 1);
184
185 if (right > median)
186 kdT_Build(heap, heapIdx, heapXdi, minRight, max, median + 1, right,
187 (root << 1) + 1);
188 }
189
190
191
192 void kdT_BuildPhotonMap (struct PhotonMap *pmap)
193 {
194 Photon *nodes;
195 unsigned long i;
196 unsigned long *heapIdx, /* Photon index array */
197 *heapXdi; /* Reverse index to heapIdx */
198
199 /* Allocate kd-tree nodes and load photons from heap file */
200 if (!(nodes = calloc(pmap -> numPhotons, sizeof(Photon))))
201 error(SYSTEM, "failed in-core heap allocation in kdT_BuildPhotonMap");
202
203 rewind(pmap -> heap);
204 if (fread(nodes, sizeof(Photon), pmap -> numPhotons, pmap -> heap) !=
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) && pmap -> srcContrib) {
352 /* Lookup in contribution photon map */
353 OBJREC *srcMod;
354 const int srcIdx = photonSrcIdx(pmap, p);
355
356 if (srcIdx < 0 || srcIdx >= nsources)
357 error(INTERNAL, "invalid light source index in photon map");
358
359 srcMod = findmaterial(source [srcIdx].so);
360
361 /* Reject photon if contributions from light source which emitted it
362 * are not sought */
363 if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
364 return;
365
366 /* Reject non-caustic photon if lookup for caustic contribs */
367 if (pmap -> lookupCaustic & !p -> caustic)
368 return;
369 }
370
371 /* Squared distance to current photon (note dist2() requires doubles) */
372 VSUB(dv, pos, p -> pos);
373 d2 = DOT(dv, dv);
374
375 /* Accept photon if closer than current max dist & add to priority queue */
376 if (d2 < pmap -> maxDist2) {
377 if (pmap -> squeue.tail < sqSize) {
378 /* Priority queue not full; append photon and restore heap */
379 i = ++pmap -> squeue.tail;
380
381 while (i > 1 && sq [(i >> 1) - 1].dist2 <= d2) {
382 sq [i - 1].idx = sq [(i >> 1) - 1].idx;
383 sq [i - 1].dist2 = sq [(i >> 1) - 1].dist2;
384 i >>= 1;
385 }
386
387 sq [--i].idx = (PhotonIdx)p;
388 sq [i].dist2 = d2;
389 /* Update maxDist if we've just filled the queue */
390 if (pmap -> squeue.tail >= pmap -> squeue.len)
391 pmap -> maxDist2 = sq [0].dist2;
392 }
393 else {
394 /* Priority queue full; replace maximum, restore heap, and
395 update maxDist */
396 i = 1;
397
398 while (i <= sqSize >> 1) {
399 j = i << 1;
400 if (j < sqSize && sq [j - 1].dist2 < sq [j].dist2)
401 j++;
402 if (d2 >= sq [j - 1].dist2)
403 break;
404 sq [i - 1].idx = sq [j - 1].idx;
405 sq [i - 1].dist2 = sq [j - 1].dist2;
406 i = j;
407 }
408
409 sq [--i].idx = (PhotonIdx)p;
410 sq [i].dist2 = d2;
411 pmap -> maxDist2 = sq [0].dist2;
412 }
413 }
414 }
415
416
417
418 void kdT_FindPhotons (struct PhotonMap *pmap, const FVECT pos,
419 const FVECT norm)
420 {
421 float p [3], n [3];
422
423 /* Photon pos & normal stored at lower precision */
424 VCOPY(p, pos);
425 VCOPY(n, norm);
426 kdT_FindNearest(pmap, p, n, 1);
427 }
428
429
430
431 static void kdT_Find1Nearest (PhotonMap *pmap, const float pos [3],
432 const float norm [3], Photon **photon,
433 unsigned long node)
434 /* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to
435 * pos with similar normal. Note that all heap and queue indices are
436 * 1-based, but accesses to the arrays are 0-based! */
437 {
438 Photon *p = (Photon*)pmap -> store.nodes + node - 1;
439 /* Signed distance to current photon's splitting plane */
440 float d = pos [p -> discr] - p -> pos [p -> discr], d2 = d * d,
441 dv [3];
442
443 /* Search subtree closer to pos first; exclude other subtree if the
444 distance to the splitting plane is greater than maxDist */
445 if (d < 0) {
446 if (node << 1 <= pmap -> numPhotons)
447 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
448
449 if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)
450 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
451 }
452 else {
453 if (node << 1 < pmap -> numPhotons)
454 kdT_Find1Nearest(pmap, pos, norm, photon, (node << 1) + 1);
455
456 if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)
457 kdT_Find1Nearest(pmap, pos, norm, photon, node << 1);
458 }
459
460 /* Squared distance to current photon */
461 VSUB(dv, pos, p -> pos);
462 d2 = DOT(dv, dv);
463
464 if (d2 < pmap -> maxDist2 &&
465 DOT(norm, p -> norm) > PMAP_NORM_TOL * 127 * frandom()) {
466 /* Closest photon so far with similar normal. We allow for tolerance
467 * to account for perturbation in the latter; note the photon normal
468 * is coded in the range [-127,127], hence we factor this in */
469 pmap -> maxDist2 = d2;
470 *photon = p;
471 }
472 }
473
474
475
476 void kdT_Find1Photon (struct PhotonMap *pmap, const FVECT pos,
477 const FVECT norm, Photon *photon)
478 {
479 float p [3], n [3];
480 Photon *pnn;
481
482 /* Photon pos & normal stored at lower precision */
483 VCOPY(p, pos);
484 VCOPY(n, norm);
485 kdT_Find1Nearest(pmap, p, n, &pnn, 1);
486 memcpy(photon, pnn, sizeof(Photon));
487 }
488
489
490
491 int kdT_GetPhoton (const struct PhotonMap *pmap, PhotonIdx idx,
492 Photon *photon)
493 {
494 memcpy(photon, idx, sizeof(Photon));
495 return 0;
496 }
497
498
499
500 Photon *kdT_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
501 {
502 return idx;
503 }
504
505
506
507 PhotonIdx kdT_FirstPhoton (const struct PhotonMap* pmap)
508 {
509 return pmap -> store.nodes;
510 }
511
512
513
514 void kdT_Delete (PhotonKdTree *kdt)
515 {
516 free(kdt -> nodes);
517 kdt -> nodes = NULL;
518 }