ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapdata.c
Revision: 2.23
Committed: Wed Apr 14 11:26:25 2021 UTC (3 years, 1 month ago) by rschregle
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, HEAD
Changes since 2.22: +21 -12 lines
Log Message:
feat(mkpmap): Added -apI option for spherical ROI

File Contents

# User Rev Content
1 greg 2.11 #ifndef lint
2 rschregle 2.23 static const char RCSid[] = "$Id: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $";
3 greg 2.11 #endif
4 rschregle 2.15
5 greg 2.1 /*
6 rschregle 2.15 =========================================================================
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 greg 2.1
14     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
15     (c) Fraunhofer Institute for Solar Energy Systems,
16 rschregle 2.23 supported by the German Research Foundation
17     (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS))
18 rschregle 2.3 (c) Lucerne University of Applied Sciences and Arts,
19 rschregle 2.23 supported by the Swiss National Science Foundation
20     (SNSF #147053, "Daylight Redirecting Components")
21 rschregle 2.15 ==========================================================================
22 greg 2.1
23 rschregle 2.23 $Id: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $
24 greg 2.1 */
25    
26    
27    
28 rschregle 2.18 #include "pmapdata.h"
29 greg 2.1 #include "pmaprand.h"
30     #include "pmapmat.h"
31     #include "otypes.h"
32     #include "source.h"
33     #include "rcontrib.h"
34 rschregle 2.7 #include "random.h"
35 greg 2.1
36    
37    
38     PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
39     NULL, NULL, NULL, NULL, NULL, NULL
40     };
41    
42    
43    
44 rschregle 2.15 /* 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 rschregle 2.19 /* Ambient include/exclude set (from ambient.c) */
52     #ifndef MAXASET
53     #define MAXASET 4095
54     #endif
55     extern OBJECT ambset [MAXASET+1];
56    
57 rschregle 2.21 /* Callback to print photon attributes acc. to user defined format */
58     int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm);
59    
60 rschregle 2.15
61    
62 greg 2.1 void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
63     /* Init photon map 'n' stuff... */
64     {
65     if (!pmap)
66     return;
67    
68 rschregle 2.15 pmap -> numPhotons = 0;
69 greg 2.1 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 rschregle 2.15 pmap -> squeue.node = NULL;
79     pmap -> squeue.len = 0;
80 greg 2.1
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 rschregle 2.15 /* 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 greg 2.1 }
105    
106    
107    
108 rschregle 2.15 void initPhotonHeap (PhotonMap *pmap)
109 greg 2.1 {
110 rschregle 2.15 int fdFlags;
111 greg 2.1
112 rschregle 2.15 if (!pmap)
113     error(INTERNAL, "undefined photon map in initPhotonHeap");
114 greg 2.1
115 rschregle 2.15 if (!pmap -> heap) {
116     /* Open heap file */
117 rschregle 2.18 mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME));
118     if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b")))
119 rschregle 2.15 error(SYSTEM, "failed opening heap file in initPhotonHeap");
120 rschregle 2.18
121 greg 2.16 #ifdef F_SETFL /* XXX is there an alternate needed for Windows? */
122 rschregle 2.15 fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
123     fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND);
124 rschregle 2.18 #endif/* ftruncate(fileno(pmap -> heap), 0); */
125 greg 2.1 }
126 rschregle 2.15 }
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 rschregle 2.18 if (!pmap -> heap || !pmap -> heapBuf) {
139     /* Silently ignore undefined heap
140     error(INTERNAL, "undefined heap in flushPhotonHeap"); */
141     return;
142     }
143 rschregle 2.15
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 greg 2.1
150 rschregle 2.15 #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 greg 2.17
160 rschregle 2.18 #if NIX
161 rschregle 2.15 if (fsync(fd))
162     error(SYSTEM, "failed fsync in flushPhotonHeap");
163 greg 2.17 #endif
164 rschregle 2.18
165 rschregle 2.15 pmap -> heapBufLen = 0;
166     }
167 greg 2.1
168 rschregle 2.15
169    
170 rschregle 2.18 #ifdef DEBUG_PMAP
171 rschregle 2.15 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 greg 2.1
180 rschregle 2.15 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 greg 2.1 }
205 rschregle 2.15 #endif
206 greg 2.1
207    
208    
209 rschregle 2.15 int newPhoton (PhotonMap* pmap, const RAY* ray)
210 greg 2.1 {
211 rschregle 2.19 unsigned i;
212 rschregle 2.15 Photon photon;
213 greg 2.1 COLOR photonFlux;
214    
215     /* Account for distribution ratio */
216     if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio)
217 rschregle 2.15 return -1;
218 greg 2.1
219     /* Don't store on sources */
220     if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
221 rschregle 2.15 return -1;
222 greg 2.1
223 rschregle 2.20 /* Ignore photon if modifier in/outside exclude/include set */
224 rschregle 2.19 if (ambincl != -1 && ray -> ro &&
225     ambincl != inset(ambset, ray -> ro -> omod))
226     return -1;
227    
228 rschregle 2.23 if (pmapNumROI && pmapROI) {
229 rschregle 2.19 unsigned inROI = 0;
230 rschregle 2.23 FVECT photonDist;
231 rschregle 2.19
232 rschregle 2.23 /* 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 rschregle 2.19 if (!inROI)
246     return -1;
247 rschregle 2.18 }
248 rschregle 2.20
249 rschregle 2.19 /* Adjust flux according to distribution ratio and ray weight */
250 rschregle 2.21 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 rschregle 2.19 scalecolor(photonFlux,
256     ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
257     : 1));
258 rschregle 2.21 #else
259     scalecolor(photonFlux,
260     1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1));
261     #endif
262 rschregle 2.19 setPhotonFlux(&photon, photonFlux);
263 rschregle 2.20
264 rschregle 2.19 /* 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 rschregle 2.15
285 rschregle 2.19 if (!pmap -> heapBuf) {
286     /* Lazily allocate heap buffa */
287 rschregle 2.18 #if NIX
288 rschregle 2.19 /* 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 rschregle 2.15 #endif
297 rschregle 2.19 if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
298     error(SYSTEM, "failed heap buffer allocation in newPhoton");
299     pmap -> heapBufLen = 0;
300     }
301 rschregle 2.15
302 rschregle 2.19 /* 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 greg 2.1
309 rschregle 2.19 pmap -> numPhotons++;
310 rschregle 2.21
311     /* Print photon attributes */
312     if (printPhoton)
313     /* Non-const kludge */
314     printPhoton((RAY*)ray, &photon, pmap);
315 greg 2.1
316 rschregle 2.15 return 0;
317 greg 2.1 }
318    
319    
320    
321 rschregle 2.15 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 rschregle 2.18 char nuHeapFname [sizeof(PMAP_TMPFNAME)];
329 rschregle 2.15 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 rschregle 2.18 if (fseek(pmap -> heap, 0, SEEK_END) < 0)
339     error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap");
340 rschregle 2.15 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 rschregle 2.18
356 rschregle 2.15 /* 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 greg 2.1 }
364    
365 rschregle 2.15 /* We REALLY don't need yet another @%&*! heap just to hold the scaled
366     * photons, but can't think of a quicker fix... */
367 rschregle 2.18 mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME));
368     if (!(nuHeap = fopen(nuHeapFname, "w+b")))
369 rschregle 2.15 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 greg 2.1
378 rschregle 2.18 while (!feof(pmap -> heap)) {
379     #ifdef DEBUG_PMAP
380     printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap));
381     #endif
382 rschregle 2.15 pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
383 rschregle 2.18 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 rschregle 2.15
391 rschregle 2.18 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 rschregle 2.20 getPhotonFlux(p, flux);
412 rschregle 2.18
413     if (photonFlux) {
414     scalecolor(flux, photonFlux [isContribPmap(pmap) ?
415     photonSrcIdx(pmap, p) : 0]);
416     setPhotonFlux(p, flux);
417 greg 2.1 }
418 rschregle 2.18
419     /* Update average photon flux; need a double here */
420     addcolor(avgFlux, flux);
421     }
422 greg 2.1
423 rschregle 2.18 /* Write modified photons to new heap */
424     fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
425 rschregle 2.15
426 rschregle 2.18 if (ferror(nuHeap))
427     error(SYSTEM, "failed postprocessing photon flux in "
428     "buildPhotonMap");
429 rschregle 2.15
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 rschregle 2.18 pmap -> heapBufSize, pmap -> heap);
452 rschregle 2.15
453 rschregle 2.18 for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
454     VSUB(d, p -> pos, CoG);
455     CoGdist += DOT(d, d);
456     }
457 rschregle 2.15 }
458    
459     pmap -> CoGdist = CoGdist /= pmap -> numPhotons;
460    
461 rschregle 2.18 /* Swap heaps, discarding unscaled photons */
462 rschregle 2.15 fclose(pmap -> heap);
463 rschregle 2.18 unlink(pmap -> heapFname);
464 rschregle 2.15 pmap -> heap = nuHeap;
465 rschregle 2.18 strcpy(pmap -> heapFname, nuHeapFname);
466 rschregle 2.15
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 rschregle 2.18 unlink(pmap -> heapFname);
477 rschregle 2.15 pmap -> heap = NULL;
478     pmap -> heapBuf = NULL;
479 greg 2.1 }
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 rschregle 2.8 /* Coefficient for adaptive maximum search radius */
494     #define PMAP_MAXDIST_COEFF 100
495    
496 greg 2.1 void findPhotons (PhotonMap* pmap, const RAY* ray)
497     {
498     int redo = 0;
499    
500 rschregle 2.15 if (!pmap -> squeue.len) {
501 greg 2.1 /* Lazy init priority queue */
502 rschregle 2.15 #ifdef PMAP_OOC
503     OOC_InitFindPhotons(pmap);
504     #else
505     kdT_InitFindPhotons(pmap);
506     #endif
507 greg 2.1 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 rschregle 2.9 /* SQUARED max search radius limit is based on avg photon distance to
516 rschregle 2.8 * centre of gravity, unless fixed by user (maxDistFix > 0) */
517 rschregle 2.15 pmap -> maxDist0 = pmap -> maxDist2Limit =
518 rschregle 2.9 maxDistFix > 0 ? maxDistFix * maxDistFix
519 rschregle 2.15 : PMAP_MAXDIST_COEFF * pmap -> squeue.len *
520     pmap -> CoGdist / pmap -> numPhotons;
521 greg 2.1 }
522    
523     do {
524 rschregle 2.15 pmap -> squeue.tail = 0;
525     pmap -> maxDist2 = pmap -> maxDist0;
526 greg 2.1
527     /* Search position is ray -> rorg for volume photons, since we have no
528     intersection point. Normals are ignored -- these are incident
529 rschregle 2.22 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 greg 2.1 if (isVolumePmap(pmap)) {
535 rschregle 2.15 #ifdef PMAP_OOC
536     OOC_FindPhotons(pmap, ray -> rorg, NULL);
537     #else
538     kdT_FindPhotons(pmap, ray -> rorg, NULL);
539     #endif
540 greg 2.1 }
541     else {
542 rschregle 2.15 #ifdef PMAP_OOC
543     OOC_FindPhotons(pmap, ray -> rop, ray -> ron);
544     #else
545     kdT_FindPhotons(pmap, ray -> rop, ray -> ron);
546     #endif
547 greg 2.1 }
548 rschregle 2.3
549 rschregle 2.15 #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 greg 2.1 /* Short lookup; too few photons found */
559 rschregle 2.15 if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) {
560 greg 2.1 /* 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 rschregle 2.8 #ifdef PMAP_LOOKUP_WARN
565 greg 2.1 sprintf(errmsg,
566     "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
567 rschregle 2.15 pmap -> squeue.tail, pmap -> squeue.len,
568     pmapName [pmap -> type],
569     ray -> rop [0], ray -> rop [1], ray -> rop [2],
570 greg 2.1 ray -> ro ? ray -> ro -> oname : "<null>");
571     error(WARNING, errmsg);
572 rschregle 2.8 #endif
573 rschregle 2.3
574 rschregle 2.8 /* Bail out after warning if maxDist is fixed */
575     if (maxDistFix > 0)
576     return;
577    
578 rschregle 2.15 if (pmap -> maxDist0 < pmap -> maxDist2Limit) {
579 greg 2.1 /* Increase max search radius if below limit & redo search */
580     pmap -> maxDist0 *= PMAP_MAXDIST_INC;
581 rschregle 2.8 #ifdef PMAP_LOOKUP_REDO
582 greg 2.1 redo = 1;
583 rschregle 2.8 #endif
584     #ifdef PMAP_LOOKUP_WARN
585 greg 2.1 sprintf(errmsg,
586     redo ? "restarting photon lookup with max radius %.1e"
587     : "max photon lookup radius adjusted to %.1e",
588 rschregle 2.15 pmap -> maxDist0);
589 greg 2.1 error(WARNING, errmsg);
590 rschregle 2.8 #endif
591 greg 2.1 }
592 rschregle 2.8 #ifdef PMAP_LOOKUP_REDO
593 greg 2.1 else {
594     sprintf(errmsg, "max photon lookup radius clamped to %.1e",
595 rschregle 2.15 pmap -> maxDist0);
596 greg 2.1 error(WARNING, errmsg);
597     }
598 rschregle 2.8 #endif
599 greg 2.1 }
600    
601     /* Reset successful lookup counter */
602     pmap -> numLookups = 0;
603 rschregle 2.3 }
604 greg 2.1 else {
605 rschregle 2.8 /* Bail out after warning if maxDist is fixed */
606     if (maxDistFix > 0)
607     return;
608    
609 greg 2.1 /* 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 rschregle 2.8
618 greg 2.1 } while (redo);
619     }
620    
621    
622    
623 rschregle 2.22 Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
624 rschregle 2.15 {
625 rschregle 2.22 /* 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 greg 2.1
634 rschregle 2.22 do {
635     pmap -> maxDist2 = maxDist2_0;
636 rschregle 2.15 #ifdef PMAP_OOC
637 rschregle 2.22 found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
638 rschregle 2.15 #else
639 rschregle 2.22 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 greg 2.1 }
658    
659    
660    
661 rschregle 2.15 void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon)
662 greg 2.1 {
663 rschregle 2.15 #ifdef PMAP_OOC
664 rschregle 2.18 if (OOC_GetPhoton(pmap, idx, photon))
665 rschregle 2.15 #else
666     if (kdT_GetPhoton(pmap, idx, photon))
667     #endif
668     error(INTERNAL, "failed photon lookup");
669 greg 2.1 }
670    
671    
672    
673 rschregle 2.15 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 greg 2.1 }
681    
682    
683    
684 rschregle 2.15 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 greg 2.1 }
692    
693    
694    
695     void deletePhotons (PhotonMap* pmap)
696     {
697 rschregle 2.15 #ifdef PMAP_OOC
698     OOC_Delete(&pmap -> store);
699     #else
700     kdT_Delete(&pmap -> store);
701     #endif
702    
703     free(pmap -> squeue.node);
704 greg 2.1 free(pmap -> biasCompHist);
705    
706 rschregle 2.15 pmap -> numPhotons = pmap -> minGather = pmap -> maxGather =
707     pmap -> squeue.len = pmap -> squeue.tail = 0;
708 greg 2.1 }