1 |
#ifndef lint |
2 |
static const char RCSid[] = "$Id: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $"; |
3 |
#endif |
4 |
|
5 |
/* |
6 |
========================================================================= |
7 |
Photon map types and interface to nearest neighbour lookups in underlying |
8 |
point cloud data structure. |
9 |
|
10 |
The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). |
11 |
This can be overriden with the PMAP_OOC compiletime switch, which enables |
12 |
an out-of-core octree (see oococt.{h,c}). |
13 |
|
14 |
Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) |
15 |
(c) Fraunhofer Institute for Solar Energy Systems, |
16 |
supported by the German Research Foundation |
17 |
(DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) |
18 |
(c) Lucerne University of Applied Sciences and Arts, |
19 |
supported by the Swiss National Science Foundation |
20 |
(SNSF #147053, "Daylight Redirecting Components") |
21 |
========================================================================== |
22 |
|
23 |
$Id: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $ |
24 |
*/ |
25 |
|
26 |
|
27 |
|
28 |
#include "pmapdata.h" |
29 |
#include "pmaprand.h" |
30 |
#include "pmapmat.h" |
31 |
#include "otypes.h" |
32 |
#include "source.h" |
33 |
#include "rcontrib.h" |
34 |
#include "random.h" |
35 |
|
36 |
|
37 |
|
38 |
PhotonMap *photonMaps [NUM_PMAP_TYPES] = { |
39 |
NULL, NULL, NULL, NULL, NULL, NULL |
40 |
}; |
41 |
|
42 |
|
43 |
|
44 |
/* Include routines to handle underlying point cloud data structure */ |
45 |
#ifdef PMAP_OOC |
46 |
#include "pmapooc.c" |
47 |
#else |
48 |
#include "pmapkdt.c" |
49 |
#endif |
50 |
|
51 |
/* Ambient include/exclude set (from ambient.c) */ |
52 |
#ifndef MAXASET |
53 |
#define MAXASET 4095 |
54 |
#endif |
55 |
extern OBJECT ambset [MAXASET+1]; |
56 |
|
57 |
/* Callback to print photon attributes acc. to user defined format */ |
58 |
int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); |
59 |
|
60 |
|
61 |
|
62 |
void initPhotonMap (PhotonMap *pmap, PhotonMapType t) |
63 |
/* Init photon map 'n' stuff... */ |
64 |
{ |
65 |
if (!pmap) |
66 |
return; |
67 |
|
68 |
pmap -> numPhotons = 0; |
69 |
pmap -> biasCompHist = NULL; |
70 |
pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; |
71 |
pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; |
72 |
pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; |
73 |
pmap -> gatherTolerance = gatherTolerance; |
74 |
pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; |
75 |
pmap -> numDensity = 0; |
76 |
pmap -> distribRatio = 1; |
77 |
pmap -> type = t; |
78 |
pmap -> squeue.node = NULL; |
79 |
pmap -> squeue.len = 0; |
80 |
|
81 |
/* Init local RNG state */ |
82 |
pmap -> randState [0] = 10243; |
83 |
pmap -> randState [1] = 39829; |
84 |
pmap -> randState [2] = 9433; |
85 |
pmapSeed(randSeed, pmap -> randState); |
86 |
|
87 |
/* Set up type-specific photon lookup callback */ |
88 |
pmap -> lookup = pmapLookup [t]; |
89 |
|
90 |
/* Mark primary photon ray as unused */ |
91 |
pmap -> lastPrimary.srcIdx = -1; |
92 |
pmap -> numPrimary = 0; |
93 |
pmap -> primaries = NULL; |
94 |
|
95 |
/* Init storage */ |
96 |
pmap -> heap = NULL; |
97 |
pmap -> heapBuf = NULL; |
98 |
pmap -> heapBufLen = 0; |
99 |
#ifdef PMAP_OOC |
100 |
OOC_Null(&pmap -> store); |
101 |
#else |
102 |
kdT_Null(&pmap -> store); |
103 |
#endif |
104 |
} |
105 |
|
106 |
|
107 |
|
108 |
void initPhotonHeap (PhotonMap *pmap) |
109 |
{ |
110 |
int fdFlags; |
111 |
|
112 |
if (!pmap) |
113 |
error(INTERNAL, "undefined photon map in initPhotonHeap"); |
114 |
|
115 |
if (!pmap -> heap) { |
116 |
/* Open heap file */ |
117 |
mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); |
118 |
if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) |
119 |
error(SYSTEM, "failed opening heap file in initPhotonHeap"); |
120 |
|
121 |
#ifdef F_SETFL /* XXX is there an alternate needed for Windows? */ |
122 |
fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); |
123 |
fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); |
124 |
#endif/* ftruncate(fileno(pmap -> heap), 0); */ |
125 |
} |
126 |
} |
127 |
|
128 |
|
129 |
|
130 |
void flushPhotonHeap (PhotonMap *pmap) |
131 |
{ |
132 |
int fd; |
133 |
const unsigned long len = pmap -> heapBufLen * sizeof(Photon); |
134 |
|
135 |
if (!pmap) |
136 |
error(INTERNAL, "undefined photon map in flushPhotonHeap"); |
137 |
|
138 |
if (!pmap -> heap || !pmap -> heapBuf) { |
139 |
/* Silently ignore undefined heap |
140 |
error(INTERNAL, "undefined heap in flushPhotonHeap"); */ |
141 |
return; |
142 |
} |
143 |
|
144 |
/* Atomically seek and write block to heap */ |
145 |
/* !!! Unbuffered I/O via pwrite() avoids potential race conditions |
146 |
* !!! and buffer corruption which can occur with lseek()/fseek() |
147 |
* !!! followed by write()/fwrite(). */ |
148 |
fd = fileno(pmap -> heap); |
149 |
|
150 |
#ifdef DEBUG_PMAP |
151 |
sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), |
152 |
pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon)); |
153 |
eputs(errmsg); |
154 |
#endif |
155 |
|
156 |
/*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ |
157 |
if (write(fd, pmap -> heapBuf, len) != len) |
158 |
error(SYSTEM, "failed append to heap file in flushPhotonHeap"); |
159 |
|
160 |
#if NIX |
161 |
if (fsync(fd)) |
162 |
error(SYSTEM, "failed fsync in flushPhotonHeap"); |
163 |
#endif |
164 |
|
165 |
pmap -> heapBufLen = 0; |
166 |
} |
167 |
|
168 |
|
169 |
|
170 |
#ifdef DEBUG_PMAP |
171 |
static int checkPhotonHeap (FILE *file) |
172 |
/* Check heap for nonsensical or duplicate photons */ |
173 |
{ |
174 |
Photon p, lastp; |
175 |
int i, dup; |
176 |
|
177 |
rewind(file); |
178 |
memset(&lastp, 0, sizeof(lastp)); |
179 |
|
180 |
while (fread(&p, sizeof(p), 1, file)) { |
181 |
dup = 1; |
182 |
|
183 |
for (i = 0; i <= 2; i++) { |
184 |
if (p.pos [i] < thescene.cuorg [i] || |
185 |
p.pos [i] > thescene.cuorg [i] + thescene.cusize) { |
186 |
|
187 |
sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", |
188 |
p.pos [0], p.pos [1], p.pos [2]); |
189 |
error(WARNING, errmsg); |
190 |
} |
191 |
|
192 |
dup &= p.pos [i] == lastp.pos [i]; |
193 |
} |
194 |
|
195 |
if (dup) { |
196 |
sprintf(errmsg, |
197 |
"consecutive duplicate photon in heap at [%f, %f, %f]\n", |
198 |
p.pos [0], p.pos [1], p.pos [2]); |
199 |
error(WARNING, errmsg); |
200 |
} |
201 |
} |
202 |
|
203 |
return 0; |
204 |
} |
205 |
#endif |
206 |
|
207 |
|
208 |
|
209 |
int newPhoton (PhotonMap* pmap, const RAY* ray) |
210 |
{ |
211 |
unsigned i; |
212 |
Photon photon; |
213 |
COLOR photonFlux; |
214 |
|
215 |
/* Account for distribution ratio */ |
216 |
if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) |
217 |
return -1; |
218 |
|
219 |
/* Don't store on sources */ |
220 |
if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) |
221 |
return -1; |
222 |
|
223 |
/* Ignore photon if modifier in/outside exclude/include set */ |
224 |
if (ambincl != -1 && ray -> ro && |
225 |
ambincl != inset(ambset, ray -> ro -> omod)) |
226 |
return -1; |
227 |
|
228 |
if (pmapNumROI && pmapROI) { |
229 |
unsigned inROI = 0; |
230 |
FVECT photonDist; |
231 |
|
232 |
/* Store photon if within a region of interest (for ze Ecksperts!) |
233 |
Note size of spherical ROI is squared. */ |
234 |
for (i = 0; !inROI && i < pmapNumROI; i++) { |
235 |
VSUB(photonDist, ray -> rop, pmapROI [i].pos); |
236 |
|
237 |
inROI = ( |
238 |
PMAP_ROI_ISSPHERE(pmapROI + i) |
239 |
? DOT(photonDist, photonDist) <= pmapROI [i].siz [0] |
240 |
: fabs(photonDist [0]) <= pmapROI [i].siz [0] && |
241 |
fabs(photonDist [1]) <= pmapROI [i].siz [1] && |
242 |
fabs(photonDist [2]) <= pmapROI [i].siz [2] |
243 |
); |
244 |
} |
245 |
if (!inROI) |
246 |
return -1; |
247 |
} |
248 |
|
249 |
/* Adjust flux according to distribution ratio and ray weight */ |
250 |
copycolor(photonFlux, ray -> rcol); |
251 |
#if 0 |
252 |
/* Factored out ray -> rweight as deprecated (?) for pmap, and infact |
253 |
erroneously attenuates volume photon flux based on extinction, |
254 |
which is already factored in by photonParticipate() */ |
255 |
scalecolor(photonFlux, |
256 |
ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio |
257 |
: 1)); |
258 |
#else |
259 |
scalecolor(photonFlux, |
260 |
1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1)); |
261 |
#endif |
262 |
setPhotonFlux(&photon, photonFlux); |
263 |
|
264 |
/* Set photon position and flags */ |
265 |
VCOPY(photon.pos, ray -> rop); |
266 |
photon.flags = 0; |
267 |
photon.caustic = PMAP_CAUSTICRAY(ray); |
268 |
|
269 |
/* Set contrib photon's primary ray and subprocess index (the latter |
270 |
* to linearise the primary ray indices after photon distribution is |
271 |
* complete). Also set primary ray's source index, thereby marking it |
272 |
* as used. */ |
273 |
if (isContribPmap(pmap)) { |
274 |
photon.primary = pmap -> numPrimary; |
275 |
photon.proc = PMAP_GETRAYPROC(ray); |
276 |
pmap -> lastPrimary.srcIdx = ray -> rsrc; |
277 |
} |
278 |
else photon.primary = 0; |
279 |
|
280 |
/* Set normal */ |
281 |
for (i = 0; i <= 2; i++) |
282 |
photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i] |
283 |
: ray -> ron [i]); |
284 |
|
285 |
if (!pmap -> heapBuf) { |
286 |
/* Lazily allocate heap buffa */ |
287 |
#if NIX |
288 |
/* Randomise buffa size to temporally decorellate flushes in |
289 |
* multiprocessing mode */ |
290 |
srandom(randSeed + getpid()); |
291 |
pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom()); |
292 |
#else |
293 |
/* Randomisation disabled for single processes on WIN; also useful |
294 |
* for reproducability during debugging */ |
295 |
pmap -> heapBufSize = PMAP_HEAPBUFSIZE; |
296 |
#endif |
297 |
if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) |
298 |
error(SYSTEM, "failed heap buffer allocation in newPhoton"); |
299 |
pmap -> heapBufLen = 0; |
300 |
} |
301 |
|
302 |
/* Photon initialised; now append to heap buffa */ |
303 |
memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); |
304 |
|
305 |
if (++pmap -> heapBufLen >= pmap -> heapBufSize) |
306 |
/* Heap buffa full, flush to heap file */ |
307 |
flushPhotonHeap(pmap); |
308 |
|
309 |
pmap -> numPhotons++; |
310 |
|
311 |
/* Print photon attributes */ |
312 |
if (printPhoton) |
313 |
/* Non-const kludge */ |
314 |
printPhoton((RAY*)ray, &photon, pmap); |
315 |
|
316 |
return 0; |
317 |
} |
318 |
|
319 |
|
320 |
|
321 |
void buildPhotonMap (PhotonMap *pmap, double *photonFlux, |
322 |
PhotonPrimaryIdx *primaryOfs, unsigned nproc) |
323 |
{ |
324 |
unsigned long n, nCheck = 0; |
325 |
unsigned i; |
326 |
Photon *p; |
327 |
COLOR flux; |
328 |
char nuHeapFname [sizeof(PMAP_TMPFNAME)]; |
329 |
FILE *nuHeap; |
330 |
/* Need double here to reduce summation errors */ |
331 |
double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0; |
332 |
FVECT d; |
333 |
|
334 |
if (!pmap) |
335 |
error(INTERNAL, "undefined photon map in buildPhotonMap"); |
336 |
|
337 |
/* Get number of photons from heapfile size */ |
338 |
if (fseek(pmap -> heap, 0, SEEK_END) < 0) |
339 |
error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap"); |
340 |
pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); |
341 |
|
342 |
if (!pmap -> numPhotons) |
343 |
error(INTERNAL, "empty photon map in buildPhotonMap"); |
344 |
|
345 |
if (!pmap -> heap) |
346 |
error(INTERNAL, "no heap in buildPhotonMap"); |
347 |
|
348 |
#ifdef DEBUG_PMAP |
349 |
eputs("Checking photon heap consistency...\n"); |
350 |
checkPhotonHeap(pmap -> heap); |
351 |
|
352 |
sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); |
353 |
eputs(errmsg); |
354 |
#endif |
355 |
|
356 |
/* Allocate heap buffa */ |
357 |
if (!pmap -> heapBuf) { |
358 |
pmap -> heapBufSize = PMAP_HEAPBUFSIZE; |
359 |
pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); |
360 |
if (!pmap -> heapBuf) |
361 |
error(SYSTEM, "failed to allocate postprocessed photon heap in" |
362 |
"buildPhotonMap"); |
363 |
} |
364 |
|
365 |
/* We REALLY don't need yet another @%&*! heap just to hold the scaled |
366 |
* photons, but can't think of a quicker fix... */ |
367 |
mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); |
368 |
if (!(nuHeap = fopen(nuHeapFname, "w+b"))) |
369 |
error(SYSTEM, "failed to open postprocessed photon heap in " |
370 |
"buildPhotonMap"); |
371 |
|
372 |
rewind(pmap -> heap); |
373 |
|
374 |
#ifdef DEBUG_PMAP |
375 |
eputs("Postprocessing photons...\n"); |
376 |
#endif |
377 |
|
378 |
while (!feof(pmap -> heap)) { |
379 |
#ifdef DEBUG_PMAP |
380 |
printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap)); |
381 |
#endif |
382 |
pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), |
383 |
pmap -> heapBufSize, pmap -> heap); |
384 |
#ifdef DEBUG_PMAP |
385 |
printf("Got %lu\n", pmap -> heapBufLen); |
386 |
#endif |
387 |
|
388 |
if (ferror(pmap -> heap)) |
389 |
error(SYSTEM, "failed to read photon heap in buildPhotonMap"); |
390 |
|
391 |
for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { |
392 |
/* Update min and max pos and set photon flux */ |
393 |
for (i = 0; i <= 2; i++) { |
394 |
if (p -> pos [i] < pmap -> minPos [i]) |
395 |
pmap -> minPos [i] = p -> pos [i]; |
396 |
else if (p -> pos [i] > pmap -> maxPos [i]) |
397 |
pmap -> maxPos [i] = p -> pos [i]; |
398 |
|
399 |
/* Update centre of gravity with photon position */ |
400 |
CoG [i] += p -> pos [i]; |
401 |
} |
402 |
|
403 |
if (primaryOfs) |
404 |
/* Linearise photon primary index from subprocess index using the |
405 |
* per-subprocess offsets in primaryOfs */ |
406 |
p -> primary += primaryOfs [p -> proc]; |
407 |
|
408 |
/* Scale photon's flux (hitherto normalised to 1 over RGB); in |
409 |
* case of a contrib photon map, this is done per light source, |
410 |
* and photonFlux is assumed to be an array */ |
411 |
getPhotonFlux(p, flux); |
412 |
|
413 |
if (photonFlux) { |
414 |
scalecolor(flux, photonFlux [isContribPmap(pmap) ? |
415 |
photonSrcIdx(pmap, p) : 0]); |
416 |
setPhotonFlux(p, flux); |
417 |
} |
418 |
|
419 |
/* Update average photon flux; need a double here */ |
420 |
addcolor(avgFlux, flux); |
421 |
} |
422 |
|
423 |
/* Write modified photons to new heap */ |
424 |
fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); |
425 |
|
426 |
if (ferror(nuHeap)) |
427 |
error(SYSTEM, "failed postprocessing photon flux in " |
428 |
"buildPhotonMap"); |
429 |
|
430 |
nCheck += pmap -> heapBufLen; |
431 |
} |
432 |
|
433 |
#ifdef DEBUG_PMAP |
434 |
if (nCheck < pmap -> numPhotons) |
435 |
error(INTERNAL, "truncated photon heap in buildPhotonMap"); |
436 |
#endif |
437 |
|
438 |
/* Finalise average photon flux */ |
439 |
scalecolor(avgFlux, 1.0 / pmap -> numPhotons); |
440 |
copycolor(pmap -> photonFlux, avgFlux); |
441 |
|
442 |
/* Average photon positions to get centre of gravity */ |
443 |
for (i = 0; i < 3; i++) |
444 |
pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; |
445 |
|
446 |
rewind(pmap -> heap); |
447 |
|
448 |
/* Compute average photon distance to centre of gravity */ |
449 |
while (!feof(pmap -> heap)) { |
450 |
pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), |
451 |
pmap -> heapBufSize, pmap -> heap); |
452 |
|
453 |
for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { |
454 |
VSUB(d, p -> pos, CoG); |
455 |
CoGdist += DOT(d, d); |
456 |
} |
457 |
} |
458 |
|
459 |
pmap -> CoGdist = CoGdist /= pmap -> numPhotons; |
460 |
|
461 |
/* Swap heaps, discarding unscaled photons */ |
462 |
fclose(pmap -> heap); |
463 |
unlink(pmap -> heapFname); |
464 |
pmap -> heap = nuHeap; |
465 |
strcpy(pmap -> heapFname, nuHeapFname); |
466 |
|
467 |
#ifdef PMAP_OOC |
468 |
OOC_BuildPhotonMap(pmap, nproc); |
469 |
#else |
470 |
kdT_BuildPhotonMap(pmap); |
471 |
#endif |
472 |
|
473 |
/* Trash heap and its buffa */ |
474 |
free(pmap -> heapBuf); |
475 |
fclose(pmap -> heap); |
476 |
unlink(pmap -> heapFname); |
477 |
pmap -> heap = NULL; |
478 |
pmap -> heapBuf = NULL; |
479 |
} |
480 |
|
481 |
|
482 |
|
483 |
/* Dynamic max photon search radius increase and reduction factors */ |
484 |
#define PMAP_MAXDIST_INC 4 |
485 |
#define PMAP_MAXDIST_DEC 0.9 |
486 |
|
487 |
/* Num successful lookups before reducing in max search radius */ |
488 |
#define PMAP_MAXDIST_CNT 1000 |
489 |
|
490 |
/* Threshold below which we assume increasing max radius won't help */ |
491 |
#define PMAP_SHORT_LOOKUP_THRESH 1 |
492 |
|
493 |
/* Coefficient for adaptive maximum search radius */ |
494 |
#define PMAP_MAXDIST_COEFF 100 |
495 |
|
496 |
void findPhotons (PhotonMap* pmap, const RAY* ray) |
497 |
{ |
498 |
int redo = 0; |
499 |
|
500 |
if (!pmap -> squeue.len) { |
501 |
/* Lazy init priority queue */ |
502 |
#ifdef PMAP_OOC |
503 |
OOC_InitFindPhotons(pmap); |
504 |
#else |
505 |
kdT_InitFindPhotons(pmap); |
506 |
#endif |
507 |
pmap -> minGathered = pmap -> maxGather; |
508 |
pmap -> maxGathered = pmap -> minGather; |
509 |
pmap -> totalGathered = 0; |
510 |
pmap -> numLookups = pmap -> numShortLookups = 0; |
511 |
pmap -> shortLookupPct = 0; |
512 |
pmap -> minError = FHUGE; |
513 |
pmap -> maxError = -FHUGE; |
514 |
pmap -> rmsError = 0; |
515 |
/* SQUARED max search radius limit is based on avg photon distance to |
516 |
* centre of gravity, unless fixed by user (maxDistFix > 0) */ |
517 |
pmap -> maxDist0 = pmap -> maxDist2Limit = |
518 |
maxDistFix > 0 ? maxDistFix * maxDistFix |
519 |
: PMAP_MAXDIST_COEFF * pmap -> squeue.len * |
520 |
pmap -> CoGdist / pmap -> numPhotons; |
521 |
} |
522 |
|
523 |
do { |
524 |
pmap -> squeue.tail = 0; |
525 |
pmap -> maxDist2 = pmap -> maxDist0; |
526 |
|
527 |
/* Search position is ray -> rorg for volume photons, since we have no |
528 |
intersection point. Normals are ignored -- these are incident |
529 |
directions). */ |
530 |
/* NOTE: status returned by XXX_FindPhotons() is currently ignored; |
531 |
if no photons are found, an empty queue is returned under the |
532 |
assumption all photons are too distant to contribute significant |
533 |
flux. */ |
534 |
if (isVolumePmap(pmap)) { |
535 |
#ifdef PMAP_OOC |
536 |
OOC_FindPhotons(pmap, ray -> rorg, NULL); |
537 |
#else |
538 |
kdT_FindPhotons(pmap, ray -> rorg, NULL); |
539 |
#endif |
540 |
} |
541 |
else { |
542 |
#ifdef PMAP_OOC |
543 |
OOC_FindPhotons(pmap, ray -> rop, ray -> ron); |
544 |
#else |
545 |
kdT_FindPhotons(pmap, ray -> rop, ray -> ron); |
546 |
#endif |
547 |
} |
548 |
|
549 |
#ifdef PMAP_LOOKUP_INFO |
550 |
fprintf(stderr, "%d/%d %s photons found within radius %.3f " |
551 |
"at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, |
552 |
pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), |
553 |
ray -> rop [0], ray -> rop [1], ray -> rop [2], |
554 |
ray -> ro ? ray -> ro -> oname : "<null>"); |
555 |
#endif |
556 |
|
557 |
if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) { |
558 |
/* Short lookup; too few photons found */ |
559 |
if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { |
560 |
/* Ignore short lookups which return fewer than |
561 |
* PMAP_SHORT_LOOKUP_THRESH photons under the assumption there |
562 |
* really are no photons in the vicinity, and increasing the max |
563 |
* search radius therefore won't help */ |
564 |
#ifdef PMAP_LOOKUP_WARN |
565 |
sprintf(errmsg, |
566 |
"%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", |
567 |
pmap -> squeue.tail, pmap -> squeue.len, |
568 |
pmapName [pmap -> type], |
569 |
ray -> rop [0], ray -> rop [1], ray -> rop [2], |
570 |
ray -> ro ? ray -> ro -> oname : "<null>"); |
571 |
error(WARNING, errmsg); |
572 |
#endif |
573 |
|
574 |
/* Bail out after warning if maxDist is fixed */ |
575 |
if (maxDistFix > 0) |
576 |
return; |
577 |
|
578 |
if (pmap -> maxDist0 < pmap -> maxDist2Limit) { |
579 |
/* Increase max search radius if below limit & redo search */ |
580 |
pmap -> maxDist0 *= PMAP_MAXDIST_INC; |
581 |
#ifdef PMAP_LOOKUP_REDO |
582 |
redo = 1; |
583 |
#endif |
584 |
#ifdef PMAP_LOOKUP_WARN |
585 |
sprintf(errmsg, |
586 |
redo ? "restarting photon lookup with max radius %.1e" |
587 |
: "max photon lookup radius adjusted to %.1e", |
588 |
pmap -> maxDist0); |
589 |
error(WARNING, errmsg); |
590 |
#endif |
591 |
} |
592 |
#ifdef PMAP_LOOKUP_REDO |
593 |
else { |
594 |
sprintf(errmsg, "max photon lookup radius clamped to %.1e", |
595 |
pmap -> maxDist0); |
596 |
error(WARNING, errmsg); |
597 |
} |
598 |
#endif |
599 |
} |
600 |
|
601 |
/* Reset successful lookup counter */ |
602 |
pmap -> numLookups = 0; |
603 |
} |
604 |
else { |
605 |
/* Bail out after warning if maxDist is fixed */ |
606 |
if (maxDistFix > 0) |
607 |
return; |
608 |
|
609 |
/* Increment successful lookup counter and reduce max search radius if |
610 |
* wraparound */ |
611 |
pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; |
612 |
if (!pmap -> numLookups) |
613 |
pmap -> maxDist0 *= PMAP_MAXDIST_DEC; |
614 |
|
615 |
redo = 0; |
616 |
} |
617 |
|
618 |
} while (redo); |
619 |
} |
620 |
|
621 |
|
622 |
|
623 |
Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon) |
624 |
{ |
625 |
/* Init (squared) search radius to avg photon dist to centre of gravity */ |
626 |
float maxDist2_0 = pmap -> CoGdist; |
627 |
int found = 0; |
628 |
#ifdef PMAP_LOOKUP_REDO |
629 |
#define REDO 1 |
630 |
#else |
631 |
#define REDO 0 |
632 |
#endif |
633 |
|
634 |
do { |
635 |
pmap -> maxDist2 = maxDist2_0; |
636 |
#ifdef PMAP_OOC |
637 |
found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon); |
638 |
#else |
639 |
found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); |
640 |
#endif |
641 |
if (found < 0) { |
642 |
/* Expand search radius to retry */ |
643 |
maxDist2_0 *= 2; |
644 |
#ifdef PMAP_LOOKUP_WARN |
645 |
sprintf(errmsg, "failed 1-NN photon lookup" |
646 |
#ifdef PMAP_LOOKUP_REDO |
647 |
", retrying with search radius %.2f", maxDist2_0 |
648 |
#endif |
649 |
); |
650 |
error(WARNING, errmsg); |
651 |
#endif |
652 |
} |
653 |
} while (REDO && found < 0); |
654 |
|
655 |
/* Return photon buffer containing valid photon, else NULL */ |
656 |
return found < 0 ? NULL : photon; |
657 |
} |
658 |
|
659 |
|
660 |
|
661 |
void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) |
662 |
{ |
663 |
#ifdef PMAP_OOC |
664 |
if (OOC_GetPhoton(pmap, idx, photon)) |
665 |
#else |
666 |
if (kdT_GetPhoton(pmap, idx, photon)) |
667 |
#endif |
668 |
error(INTERNAL, "failed photon lookup"); |
669 |
} |
670 |
|
671 |
|
672 |
|
673 |
Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) |
674 |
{ |
675 |
#ifdef PMAP_OOC |
676 |
return OOC_GetNearestPhoton(squeue, idx); |
677 |
#else |
678 |
return kdT_GetNearestPhoton(squeue, idx); |
679 |
#endif |
680 |
} |
681 |
|
682 |
|
683 |
|
684 |
PhotonIdx firstPhoton (const PhotonMap *pmap) |
685 |
{ |
686 |
#ifdef PMAP_OOC |
687 |
return OOC_FirstPhoton(pmap); |
688 |
#else |
689 |
return kdT_FirstPhoton(pmap); |
690 |
#endif |
691 |
} |
692 |
|
693 |
|
694 |
|
695 |
void deletePhotons (PhotonMap* pmap) |
696 |
{ |
697 |
#ifdef PMAP_OOC |
698 |
OOC_Delete(&pmap -> store); |
699 |
#else |
700 |
kdT_Delete(&pmap -> store); |
701 |
#endif |
702 |
|
703 |
free(pmap -> squeue.node); |
704 |
free(pmap -> biasCompHist); |
705 |
|
706 |
pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = |
707 |
pmap -> squeue.len = pmap -> squeue.tail = 0; |
708 |
} |