ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapdata.c
(Generate patch)

Comparing ray/src/rt/pmapdata.c (file contents):
Revision 2.3 by rschregle, Fri May 8 13:20:23 2015 UTC vs.
Revision 2.17 by greg, Thu Sep 29 20:51:07 2016 UTC

# Line 1 | Line 1
1 + #ifndef lint
2 + static const char RCSid[] = "$Id$";
3 + #endif
4 +
5   /*
6 <   ==================================================================
7 <   Photon map data structures and kd-tree handling
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     (c) Lucerne University of Applied Sciences and Arts,
17 <   supported by the Swiss National Science Foundation (SNSF, #147053)
18 <   ==================================================================  
17 >       supported by the Swiss National Science Foundation (SNSF, #147053)
18 >   ==========================================================================
19    
20     $Id$
21   */
# Line 19 | Line 28
28   #include "otypes.h"
29   #include "source.h"
30   #include "rcontrib.h"
31 + #include "random.h"
32  
33  
34  
# Line 28 | Line 38 | PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
38  
39  
40  
41 + /* Include routines to handle underlying point cloud data structure */
42 + #ifdef PMAP_OOC
43 +   #include "pmapooc.c"
44 + #else
45 +   #include "pmapkdt.c"
46 + #endif
47 +
48 +
49 +
50   void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
51   /* Init photon map 'n' stuff... */
52   {
53     if (!pmap)
54        return;
55        
56 <   pmap -> heapSize = pmap -> heapEnd = 0;
38 <   pmap -> heap = NULL;
39 <   pmap -> squeue = NULL;
56 >   pmap -> numPhotons = 0;
57     pmap -> biasCompHist = NULL;
58     pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE;
59     pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE;
# Line 46 | Line 63 | void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
63     pmap -> numDensity = 0;
64     pmap -> distribRatio = 1;
65     pmap -> type = t;
66 +   pmap -> squeue.node = NULL;
67 +   pmap -> squeue.len = 0;  
68  
69     /* Init local RNG state */
70     pmap -> randState [0] = 10243;
# Line 57 | Line 76 | void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
76     /* Set up type-specific photon lookup callback */
77     pmap -> lookup = pmapLookup [t];
78  
79 <   pmap -> primary = NULL;
80 <   pmap -> primarySize = pmap -> primaryEnd = 0;
79 >   /* Mark primary photon ray as unused */
80 >   pmap -> lastPrimary.srcIdx = -1;
81 >   pmap -> numPrimary = 0;  
82 >   pmap -> primaries = NULL;
83 >  
84 >   /* Init storage */
85 >   pmap -> heap = NULL;
86 >   pmap -> heapBuf = NULL;
87 >   pmap -> heapBufLen = 0;
88 > #ifdef PMAP_OOC
89 >   OOC_Null(&pmap -> store);
90 > #else
91 >   kdT_Null(&pmap -> store);
92 > #endif
93   }
94  
95  
96  
97 < const PhotonPrimary* addPhotonPrimary (PhotonMap *pmap, const RAY *ray)
97 > void initPhotonHeap (PhotonMap *pmap)
98   {
99 <   PhotonPrimary *prim = NULL;
99 >   int fdFlags;
100    
101 <   if (!pmap || !ray)
102 <      return NULL;
101 >   if (!pmap)
102 >      error(INTERNAL, "undefined photon map in initPhotonHeap");
103        
104 <   /* Check if last primary ray has spawned photons (srcIdx >= 0, see
105 <    * addPhoton()), in which case we keep it and allocate a new one;
106 <    * otherwise we overwrite the unused entry */
107 <   if (pmap -> primary && pmap -> primary [pmap -> primaryEnd].srcIdx >= 0)
108 <      pmap -> primaryEnd++;
109 <      
110 <   if (!pmap -> primarySize || pmap -> primaryEnd >= pmap -> primarySize) {
111 <      /* Allocate/enlarge array */
112 <      pmap -> primarySize += pmap -> heapSizeInc;
82 <      
83 <      /* Counter wraparound? */
84 <      if (pmap -> primarySize < pmap -> heapSizeInc)
85 <         error(INTERNAL, "photon primary overflow");
86 <        
87 <      pmap -> primary = (PhotonPrimary *)realloc(pmap -> primary,
88 <                                                 sizeof(PhotonPrimary) *
89 <                                                 pmap -> primarySize);
90 <      if (!pmap -> primary)
91 <         error(USER, "can't allocate photon primaries");
104 >   if (!pmap -> heap) {
105 >      /* Open heap file */
106 >      if (!(pmap -> heap = tmpfile()))
107 >         error(SYSTEM, "failed opening heap file in initPhotonHeap");
108 > #ifdef F_SETFL  /* XXX is there an alternate needed for Windows? */
109 >      fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
110 >      fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND);
111 > #endif
112 > /*      ftruncate(fileno(pmap -> heap), 0); */
113     }
114 + }
115 +
116 +
117 +
118 + void flushPhotonHeap (PhotonMap *pmap)
119 + {
120 +   int                  fd;
121 +   const unsigned long  len = pmap -> heapBufLen * sizeof(Photon);
122    
123 <   prim = pmap -> primary + pmap -> primaryEnd;
123 >   if (!pmap)
124 >      error(INTERNAL, "undefined photon map in flushPhotonHeap");
125 >
126 >   if (!pmap -> heap || !pmap -> heapBuf)
127 >      error(INTERNAL, "undefined heap in flushPhotonHeap");
128 >
129 >   /* Atomically seek and write block to heap */
130 >   /* !!! Unbuffered I/O via pwrite() avoids potential race conditions
131 >    * !!! and buffer corruption which can occur with lseek()/fseek()
132 >    * !!! followed by write()/fwrite(). */  
133 >   fd = fileno(pmap -> heap);
134    
135 <   /* Mark unused with negative source index until path spawns a photon (see
136 <    * addPhoton()) */
137 <   prim -> srcIdx = -1;
138 <      
139 <   /* Reverse incident direction to point to light source */
101 <   prim -> dir [0] = -ray -> rdir [0];
102 <   prim -> dir [1] = -ray -> rdir [1];
103 <   prim -> dir [2] = -ray -> rdir [2];
135 > #ifdef DEBUG_PMAP
136 >   sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(),
137 >           pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon));
138 >   eputs(errmsg);
139 > #endif
140  
141 <   VCOPY(prim -> org, ray -> rorg);
141 >   /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */
142 >   if (write(fd, pmap -> heapBuf, len) != len)
143 >      error(SYSTEM, "failed append to heap file in flushPhotonHeap");
144 >
145 > #if !defined(_WIN32) && !defined(_WIN64)
146 >   if (fsync(fd))
147 >      error(SYSTEM, "failed fsync in flushPhotonHeap");
148 > #endif
149 >   pmap -> heapBufLen = 0;
150 > }
151 >
152 >
153 >
154 > #ifdef DEBUG_OOC
155 > static int checkPhotonHeap (FILE *file)
156 > /* Check heap for nonsensical or duplicate photons */
157 > {
158 >   Photon   p, lastp;
159 >   int      i, dup;
160    
161 <   return prim;
161 >   rewind(file);
162 >   memset(&lastp, 0, sizeof(lastp));
163 >  
164 >   while (fread(&p, sizeof(p), 1, file)) {
165 >      dup = 1;
166 >      
167 >      for (i = 0; i <= 2; i++) {
168 >         if (p.pos [i] < thescene.cuorg [i] ||
169 >             p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
170 >            
171 >            sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
172 >                    p.pos [0], p.pos [1], p.pos [2]);
173 >            error(WARNING, errmsg);
174 >         }
175 >        
176 >         dup &= p.pos [i] == lastp.pos [i];
177 >      }
178 >      
179 >      if (dup) {
180 >         sprintf(errmsg,
181 >                 "consecutive duplicate photon in heap at [%f, %f, %f]\n",
182 >                 p.pos [0], p.pos [1], p.pos [2]);
183 >         error(WARNING, errmsg);
184 >      }
185 >   }
186 >
187 >   return 0;
188   }
189 + #endif
190  
191  
192  
193 < const Photon* addPhoton (PhotonMap* pmap, const RAY* ray)
193 > int newPhoton (PhotonMap* pmap, const RAY* ray)
194   {
195     unsigned i;
196 <   Photon* photon = NULL;
196 >   Photon photon;
197     COLOR photonFlux;
198    
199     /* Account for distribution ratio */
200     if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio)
201 <      return NULL;
201 >      return -1;
202        
203     /* Don't store on sources */
204     if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
205 <      return NULL;
205 >      return -1;
206  
126 – #if 0
127 –   if (inContribPmap(pmap)) {
128 –      /* Adding contribution photon */
129 –      if (ray -> parent && ray -> parent -> rtype & PRIMARY)
130 –         /* Add primary photon ray if parent is primary; by putting this
131 –          * here and checking the ray's immediate parent, we only add
132 –          * primaries that actually contribute photons, and we only add them
133 –          * once.  */
134 –         addPhotonPrimary(pmap, ray -> parent);
135 –      
136 –      /* Save index to primary ray (remains unchanged if primary already in
137 –       * array) */
138 –      primary = pmap -> primaryEnd;
139 –   }
140 – #endif
141 –  
207   #ifdef PMAP_ROI
208 <   /* Store photon if within region of interest -- for egg-spurtz only! */
208 >   /* Store photon if within region of interest -- for Ze Eckspertz only! */
209     if (ray -> rop [0] >= pmapROI [0] && ray -> rop [0] <= pmapROI [1] &&
210         ray -> rop [1] >= pmapROI [2] && ray -> rop [1] <= pmapROI [3] &&
211         ray -> rop [2] >= pmapROI [4] && ray -> rop [2] <= pmapROI [5])
212   #endif
213     {      
149 –      if (pmap -> heapEnd >= pmap -> heapSize) {
150 –         /* Enlarge heap */
151 –         pmap -> heapSize += pmap -> heapSizeInc;
152 –        
153 –         /* Counter wraparound? */
154 –         if (pmap -> heapSize < pmap -> heapSizeInc)
155 –            error(INTERNAL, "photon heap overflow");
156 –        
157 –         pmap -> heap = (Photon *)realloc(pmap -> heap,
158 –                                          sizeof(Photon) * pmap -> heapSize);
159 –         if (!pmap -> heap)
160 –            error(USER, "can't allocate photon heap");
161 –      }
162 –
163 –      photon = pmap -> heap + pmap -> heapEnd++;
164 –
214        /* Adjust flux according to distribution ratio and ray weight */
215        copycolor(photonFlux, ray -> rcol);  
216        scalecolor(photonFlux,
217                   ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
218                                                          : 1));
219 <      setPhotonFlux(photon, photonFlux);
219 >      setPhotonFlux(&photon, photonFlux);
220                
221        /* Set photon position and flags */
222 <      VCOPY(photon -> pos, ray -> rop);
223 <      photon -> flags = PMAP_CAUSTICRAY(ray) ? PMAP_CAUST_BIT : 0;
222 >      VCOPY(photon.pos, ray -> rop);
223 >      photon.flags = 0;
224 >      photon.caustic = PMAP_CAUSTICRAY(ray);
225  
226 <      /* Set primary ray index and mark as used for contrib photons */
226 >      /* Set contrib photon's primary ray and subprocess index (the latter
227 >       * to linearise the primary ray indices after photon distribution is
228 >       * complete). Also set primary ray's source index, thereby marking it
229 >       * as used. */
230        if (isContribPmap(pmap)) {
231 <         photon -> primary = pmap -> primaryEnd;
232 <         pmap -> primary [pmap -> primaryEnd].srcIdx = ray -> rsrc;
231 >         photon.primary = pmap -> numPrimary;
232 >         photon.proc = PMAP_GETRAYPROC(ray);
233 >         pmap -> lastPrimary.srcIdx = ray -> rsrc;
234        }
235 <      else photon -> primary = 0;
235 >      else photon.primary = 0;
236        
237 <      /* Update min and max positions & set normal */
238 <      for (i = 0; i <= 2; i++) {
239 <         if (photon -> pos [i] < pmap -> minPos [i])
240 <            pmap -> minPos [i] = photon -> pos [i];
241 <         if (photon -> pos [i] > pmap -> maxPos [i])
242 <            pmap -> maxPos [i] = photon -> pos [i];
243 <         photon -> norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
244 <                                                          : ray -> ron [i]);
237 >      /* Set normal */
238 >      for (i = 0; i <= 2; i++)
239 >         photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
240 >                                                       : ray -> ron [i]);
241 >
242 >      if (!pmap -> heapBuf) {
243 >         /* Lazily allocate heap buffa */
244 > #if 1        
245 >         /* Randomise buffa size to temporally decorellate buffa flushes */        
246 >         srandom(randSeed + getpid());
247 >         pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
248 > #else
249 >         /* Randomisation disabled for reproducability during debugging */        
250 >         pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
251 > #endif        
252 >         if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
253 >            error(SYSTEM, "failed heap buffer allocation in newPhoton");
254 >         pmap -> heapBufLen = 0;      
255        }
256 +
257 +      /* Photon initialised; now append to heap buffa */
258 +      memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
259 +                  
260 +      if (++pmap -> heapBufLen >= pmap -> heapBufSize)
261 +         /* Heap buffa full, flush to heap file */
262 +         flushPhotonHeap(pmap);
263 +
264 +      pmap -> numPhotons++;
265     }
266              
267 <   return photon;
267 >   return 0;
268   }
269  
270  
271  
272 < static void nearestNeighbours (PhotonMap* pmap, const float pos [3],
273 <                               const float norm [3], unsigned long node)
201 < /* Recursive part of findPhotons(..).
202 <   Note that all heap and priority queue index handling is 1-based, but
203 <   accesses to the arrays are 0-based! */
272 > void buildPhotonMap (PhotonMap *pmap, double *photonFlux,
273 >                     PhotonPrimaryIdx *primaryOfs, unsigned nproc)
274   {
275 <   Photon* p = &pmap -> heap [node - 1];
276 <   unsigned i, j;
277 <   /* Signed distance to current photon's splitting plane */
278 <   float d = pos [photonDiscr(*p)] - p -> pos [photonDiscr(*p)],
279 <         d2 = d * d;
280 <   PhotonSQNode* sq = pmap -> squeue;
281 <   const unsigned sqSize = pmap -> squeueSize;
282 <   float dv [3];
275 >   unsigned long  n, nCheck = 0;
276 >   unsigned       i;
277 >   Photon         *p;
278 >   COLOR          flux;
279 >   FILE           *nuHeap;
280 >   /* Need double here to reduce summation errors */
281 >   double         avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
282 >   FVECT          d;
283    
284 <   /* Search subtree closer to pos first; exclude other subtree if the
285 <      distance to the splitting plane is greater than maxDist */
216 <   if (d < 0) {
217 <      if (node << 1 <= pmap -> heapSize)
218 <         nearestNeighbours(pmap, pos, norm, node << 1);
219 <      if (d2 < pmap -> maxDist && node << 1 < pmap -> heapSize)
220 <         nearestNeighbours(pmap, pos, norm, (node << 1) + 1);
221 <   }
222 <   else {
223 <      if (node << 1 < pmap -> heapSize)
224 <         nearestNeighbours(pmap, pos, norm, (node << 1) + 1);
225 <      if (d2 < pmap -> maxDist && node << 1 <= pmap -> heapSize)
226 <         nearestNeighbours(pmap, pos, norm, node << 1);
227 <   }
228 <
229 <   /* Reject photon if normal faces away (ignored for volume photons) */
230 <   if (norm && DOT(norm, p -> norm) <= 0)
231 <      return;
284 >   if (!pmap)
285 >      error(INTERNAL, "undefined photon map in buildPhotonMap");
286        
287 <   if (isContribPmap(pmap) && pmap -> srcContrib) {
288 <      /* Lookup in contribution photon map */
289 <      OBJREC *srcMod;
290 <      const int srcIdx = photonSrcIdx(pmap, p);
291 <      
292 <      if (srcIdx < 0 || srcIdx >= nsources)
239 <         error(INTERNAL, "invalid light source index in photon map");
240 <      
241 <      srcMod = objptr(source [srcIdx].so -> omod);
287 >   /* Get number of photons from heapfile size */
288 >   fseek(pmap -> heap, 0, SEEK_END);      
289 >   pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);
290 >  
291 >   if (!pmap -> numPhotons)
292 >      error(INTERNAL, "empty photon map in buildPhotonMap");  
293  
294 <      /* Reject photon if contributions from light source which emitted it
295 <       * are not sought */
245 <      if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
246 <         return;
294 >   if (!pmap -> heap)
295 >      error(INTERNAL, "no heap in buildPhotonMap");
296  
297 <      /* Reject non-caustic photon if lookup for caustic contribs */
298 <      if (pmap -> lookupFlags & PMAP_CAUST_BIT & ~p -> flags)
299 <         return;
297 > #ifdef DEBUG_PMAP
298 >   eputs("Checking photon heap consistency...\n");
299 >   checkPhotonHeap(pmap -> heap);
300 >  
301 >   sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons);
302 >   eputs(errmsg);
303 > #endif
304 >    
305 >   /* Allocate heap buffa */
306 >   if (!pmap -> heapBuf) {
307 >      pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
308 >      pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon));
309 >      if (!pmap -> heapBuf)
310 >         error(SYSTEM, "failed to allocate postprocessed photon heap in"
311 >               "buildPhotonMap");
312     }
313 +
314 +   /* We REALLY don't need yet another @%&*! heap just to hold the scaled
315 +    * photons, but can't think of a quicker fix... */
316 +   if (!(nuHeap = tmpfile()))
317 +      error(SYSTEM, "failed to open postprocessed photon heap in "
318 +            "buildPhotonMap");
319 +            
320 +   rewind(pmap -> heap);
321 +
322 + #ifdef DEBUG_PMAP
323 +   eputs("Postprocessing photons...\n");
324 + #endif
325    
326 <   /* Squared distance to current photon */
327 <   dv [0] = pos [0] - p -> pos [0];
328 <   dv [1] = pos [1] - p -> pos [1];
329 <   dv [2] = pos [2] - p -> pos [2];
330 <   d2 = DOT(dv, dv);
326 >   while (!feof(pmap -> heap)) {
327 >      pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
328 >                                 PMAP_HEAPBUFSIZE, pmap -> heap);
329 >      
330 >      if (pmap -> heapBufLen) {
331 >         for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
332 >            /* Update min and max pos and set photon flux */
333 >            for (i = 0; i <= 2; i++) {
334 >               if (p -> pos [i] < pmap -> minPos [i])
335 >                  pmap -> minPos [i] = p -> pos [i];
336 >               else if (p -> pos [i] > pmap -> maxPos [i])
337 >                  pmap -> maxPos [i] = p -> pos [i];  
338  
339 <   /* Accept photon if closer than current max dist & add to priority queue */
340 <   if (d2 < pmap -> maxDist) {
341 <      if (pmap -> squeueEnd < sqSize) {
342 <         /* Priority queue not full; append photon and restore heap */
343 <         i = ++pmap -> squeueEnd;
344 <        
345 <         while (i > 1 && sq [(i >> 1) - 1].dist <= d2) {
346 <            sq [i - 1].photon = sq [(i >> 1) - 1].photon;
347 <            sq [i - 1].dist = sq [(i >> 1) - 1].dist;
348 <            i >>= 1;
339 >               /* Update centre of gravity with photon position */                
340 >               CoG [i] += p -> pos [i];                  
341 >            }  
342 >            
343 >            if (primaryOfs)
344 >               /* Linearise photon primary index from subprocess index using the
345 >                * per-subprocess offsets in primaryOfs */
346 >               p -> primary += primaryOfs [p -> proc];
347 >            
348 >            /* Scale photon's flux (hitherto normalised to 1 over RGB); in
349 >             * case of a contrib photon map, this is done per light source,
350 >             * and photonFlux is assumed to be an array */
351 >            getPhotonFlux(p, flux);            
352 >
353 >            if (photonFlux) {
354 >               scalecolor(flux, photonFlux [isContribPmap(pmap) ?
355 >                                               photonSrcIdx(pmap, p) : 0]);
356 >               setPhotonFlux(p, flux);
357 >            }
358 >
359 >            /* Update average photon flux; need a double here */
360 >            addcolor(avgFlux, flux);
361           }
362          
363 <         sq [--i].photon = p;
364 <         sq [i].dist = d2;
365 <         /* Update maxDist if we've just filled the queue */
366 <         if (pmap -> squeueEnd >= pmap -> squeueSize)
367 <            pmap -> maxDist = sq [0].dist;
363 >         /* Write modified photons to new heap */
364 >         fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
365 >                
366 >         if (ferror(nuHeap))
367 >            error(SYSTEM, "failed postprocessing photon flux in "
368 >                  "buildPhotonMap");
369        }
370 <      else {
371 <         /* Priority queue full; replace maximum, restore heap, and
279 <            update maxDist */
280 <         i = 1;
281 <        
282 <         while (i <= sqSize >> 1) {
283 <            j = i << 1;
284 <            if (j < sqSize && sq [j - 1].dist < sq [j].dist)
285 <               j++;
286 <            if (d2 >= sq [j - 1].dist)
287 <               break;
288 <            sq [i - 1].photon = sq [j - 1].photon;
289 <            sq [i - 1].dist = sq [j - 1].dist;
290 <            i = j;
291 <         }
292 <        
293 <         sq [--i].photon = p;
294 <         sq [i].dist = d2;
295 <         pmap -> maxDist = sq [0].dist;
296 <      }
370 >      
371 >      nCheck += pmap -> heapBufLen;
372     }
373 +
374 + #ifdef DEBUG_PMAP
375 +   if (nCheck < pmap -> numPhotons)
376 +      error(INTERNAL, "truncated photon heap in buildPhotonMap");
377 + #endif
378 +  
379 +   /* Finalise average photon flux */
380 +   scalecolor(avgFlux, 1.0 / pmap -> numPhotons);
381 +   copycolor(pmap -> photonFlux, avgFlux);
382 +
383 +   /* Average photon positions to get centre of gravity */
384 +   for (i = 0; i < 3; i++)
385 +      pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons;
386 +      
387 +   rewind(pmap -> heap);
388 +  
389 +   /* Compute average photon distance to centre of gravity */
390 +   while (!feof(pmap -> heap)) {
391 +      pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
392 +                                 PMAP_HEAPBUFSIZE, pmap -> heap);
393 +      
394 +      if (pmap -> heapBufLen)
395 +         for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
396 +            VSUB(d, p -> pos, CoG);
397 +            CoGdist += DOT(d, d);
398 +         }
399 +   }  
400 +
401 +   pmap -> CoGdist = CoGdist /= pmap -> numPhotons;
402 +
403 +   /* Swap heaps */
404 +   fclose(pmap -> heap);
405 +   pmap -> heap = nuHeap;
406 +  
407 + #ifdef PMAP_OOC
408 +   OOC_BuildPhotonMap(pmap, nproc);
409 + #else
410 +   /* kd-tree not parallelised */
411 +   kdT_BuildPhotonMap(pmap);
412 + #endif
413 +
414 +   /* Trash heap and its buffa */
415 +   free(pmap -> heapBuf);
416 +   fclose(pmap -> heap);
417 +   pmap -> heap = NULL;
418 +   pmap -> heapBuf = NULL;
419   }
420  
421  
# Line 309 | Line 430 | static void nearestNeighbours (PhotonMap* pmap, const
430   /* Threshold below which we assume increasing max radius won't help */
431   #define PMAP_SHORT_LOOKUP_THRESH 1
432  
433 + /* Coefficient for adaptive maximum search radius */
434 + #define PMAP_MAXDIST_COEFF 100
435 +
436   void findPhotons (PhotonMap* pmap, const RAY* ray)
437   {
314 –   float pos [3], norm [3];
438     int redo = 0;
439    
440 <   if (!pmap -> squeue) {
440 >   if (!pmap -> squeue.len) {
441        /* Lazy init priority queue */
442 <      pmap -> squeueSize = pmap -> maxGather + 1;
443 <      pmap -> squeue = (PhotonSQNode*)malloc(pmap -> squeueSize *
444 <                                             sizeof(PhotonSQNode));
445 <      if (!pmap -> squeue)
446 <         error(USER, "can't allocate photon priority queue");
324 <        
442 > #ifdef PMAP_OOC
443 >      OOC_InitFindPhotons(pmap);
444 > #else
445 >      kdT_InitFindPhotons(pmap);      
446 > #endif      
447        pmap -> minGathered = pmap -> maxGather;
448        pmap -> maxGathered = pmap -> minGather;
449        pmap -> totalGathered = 0;
# Line 330 | Line 452 | void findPhotons (PhotonMap* pmap, const RAY* ray)
452        pmap -> minError = FHUGE;
453        pmap -> maxError = -FHUGE;
454        pmap -> rmsError = 0;
455 < #ifdef PMAP_MAXDIST_ABS
456 <      /* Treat maxDistCoeff as an *absolute* and *fixed* max search radius.
457 <         Primarily intended for debugging; FOR ZE ECKSPURTZ ONLY! */
458 <      pmap -> maxDist0 = pmap -> maxDistLimit = maxDistCoeff;
459 < #else
460 <      /* Maximum search radius limit based on avg photon distance to
339 <       * centre of gravity */
340 <      pmap -> maxDist0 = pmap -> maxDistLimit =
341 <         maxDistCoeff * pmap -> squeueSize * pmap -> CoGdist /
342 <         pmap -> heapSize;
343 < #endif      
455 >      /* SQUARED max search radius limit is based on avg photon distance to
456 >       * centre of gravity, unless fixed by user (maxDistFix > 0) */
457 >      pmap -> maxDist0 = pmap -> maxDist2Limit =
458 >         maxDistFix > 0 ? maxDistFix * maxDistFix
459 >                        : PMAP_MAXDIST_COEFF * pmap -> squeue.len *
460 >                          pmap -> CoGdist / pmap -> numPhotons;
461     }
462  
463     do {
464 <      pmap -> squeueEnd = 0;
465 <      pmap -> maxDist = pmap -> maxDist0;
464 >      pmap -> squeue.tail = 0;
465 >      pmap -> maxDist2 = pmap -> maxDist0;
466          
467        /* Search position is ray -> rorg for volume photons, since we have no
468           intersection point. Normals are ignored -- these are incident
469           directions). */  
470        if (isVolumePmap(pmap)) {
471 <         VCOPY(pos, ray -> rorg);
472 <         nearestNeighbours(pmap, pos, NULL, 1);
471 > #ifdef PMAP_OOC
472 >         OOC_FindPhotons(pmap, ray -> rorg, NULL);
473 > #else
474 >         kdT_FindPhotons(pmap, ray -> rorg, NULL);
475 > #endif                  
476        }
477        else {
478 <         VCOPY(pos, ray -> rop);
479 <         VCOPY(norm, ray -> ron);
480 <         nearestNeighbours(pmap, pos, norm, 1);
478 > #ifdef PMAP_OOC
479 >         OOC_FindPhotons(pmap, ray -> rop, ray -> ron);        
480 > #else
481 >         kdT_FindPhotons(pmap, ray -> rop, ray -> ron);
482 > #endif
483        }
484  
485 < #ifndef PMAP_MAXDIST_ABS      
486 <      if (pmap -> squeueEnd < pmap -> squeueSize * pmap -> gatherTolerance) {
485 > #ifdef PMAP_LOOKUP_INFO
486 >      fprintf(stderr, "%d/%d %s photons found within radius %.3f "
487 >              "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail,
488 >              pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2),
489 >              ray -> rop [0], ray -> rop [1], ray -> rop [2],
490 >              ray -> ro ? ray -> ro -> oname : "<null>");
491 > #endif      
492 >            
493 >      if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) {
494           /* Short lookup; too few photons found */
495 <         if (pmap -> squeueEnd > PMAP_SHORT_LOOKUP_THRESH) {
495 >         if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) {
496              /* Ignore short lookups which return fewer than
497               * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there
498               * really are no photons in the vicinity, and increasing the max
499               * search radius therefore won't help */
500 <   #ifdef PMAP_LOOKUP_WARN
500 > #ifdef PMAP_LOOKUP_WARN
501              sprintf(errmsg,
502                      "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
503 <                    pmap -> squeueEnd, pmap -> squeueSize,
504 <                    pmapName [pmap -> type], pos [0], pos [1], pos [2],
503 >                    pmap -> squeue.tail, pmap -> squeue.len,
504 >                    pmapName [pmap -> type],
505 >                    ray -> rop [0], ray -> rop [1], ray -> rop [2],
506                      ray -> ro ? ray -> ro -> oname : "<null>");
507              error(WARNING, errmsg);        
508 <   #endif            
508 > #endif            
509  
510 <            if (pmap -> maxDist0 < pmap -> maxDistLimit) {
510 >            /* Bail out after warning if maxDist is fixed */
511 >            if (maxDistFix > 0)
512 >               return;
513 >            
514 >            if (pmap -> maxDist0 < pmap -> maxDist2Limit) {
515                 /* Increase max search radius if below limit & redo search */
516                 pmap -> maxDist0 *= PMAP_MAXDIST_INC;
517 <   #ifdef PMAP_LOOKUP_REDO
517 > #ifdef PMAP_LOOKUP_REDO
518                 redo = 1;
519 <   #endif              
520 <
387 <   #ifdef PMAP_LOOKUP_WARN
519 > #endif              
520 > #ifdef PMAP_LOOKUP_WARN
521                 sprintf(errmsg,
522                         redo ? "restarting photon lookup with max radius %.1e"
523                              : "max photon lookup radius adjusted to %.1e",
524                         pmap -> maxDist0);
525                 error(WARNING, errmsg);
526 <   #endif
526 > #endif
527              }
528 <   #ifdef PMAP_LOOKUP_REDO
528 > #ifdef PMAP_LOOKUP_REDO
529              else {
530                 sprintf(errmsg, "max photon lookup radius clamped to %.1e",
531                         pmap -> maxDist0);
532                 error(WARNING, errmsg);
533              }
534 <   #endif
534 > #endif
535           }
536          
537           /* Reset successful lookup counter */
538           pmap -> numLookups = 0;
539        }
540        else {
541 +         /* Bail out after warning if maxDist is fixed */
542 +         if (maxDistFix > 0)
543 +            return;
544 +            
545           /* Increment successful lookup counter and reduce max search radius if
546            * wraparound */
547           pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT;
# Line 413 | Line 550 | void findPhotons (PhotonMap* pmap, const RAY* ray)
550              
551           redo = 0;
552        }
553 < #endif
553 >      
554     } while (redo);
555   }
556  
557  
558  
559 < static void nearest1Neighbour (PhotonMap *pmap, const float pos [3],
423 <                               const float norm [3], Photon **photon,
424 <                               unsigned long node)
425 < /* Recursive part of find1Photon(..).
426 <   Note that all heap index handling is 1-based, but accesses to the
427 <   arrays are 0-based! */
559 > void find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
560   {
561 <   Photon *p = pmap -> heap + node - 1;
430 <   /* Signed distance to current photon's splitting plane */
431 <   float d = pos [photonDiscr(*p)] - p -> pos [photonDiscr(*p)],
432 <         d2 = d * d;
433 <   float dv [3];
561 >   pmap -> maxDist2 = thescene.cusize;  /* ? */
562    
563 <   /* Search subtree closer to pos first; exclude other subtree if the
564 <      distance to the splitting plane is greater than maxDist */
565 <   if (d < 0) {
566 <      if (node << 1 <= pmap -> heapSize)
567 <         nearest1Neighbour(pmap, pos, norm, photon, node << 1);
440 <      if (d2 < pmap -> maxDist && node << 1 < pmap -> heapSize)
441 <         nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
442 <   }
443 <   else {
444 <      if (node << 1 < pmap -> heapSize)
445 <         nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
446 <      if (d2 < pmap -> maxDist && node << 1 <= pmap -> heapSize)
447 <         nearest1Neighbour(pmap, pos, norm, photon, node << 1);
448 <   }
449 <  
450 <   /* Squared distance to current photon */
451 <   dv [0] = pos [0] - p -> pos [0];
452 <   dv [1] = pos [1] - p -> pos [1];
453 <   dv [2] = pos [2] - p -> pos [2];
454 <   d2 = DOT(dv, dv);
455 <  
456 <   if (d2 < pmap -> maxDist && DOT(norm, p -> norm) > 0) {
457 <      /* Closest photon so far with similar normal */
458 <      pmap -> maxDist = d2;
459 <      *photon = p;
460 <   }
563 > #ifdef PMAP_OOC
564 >   OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
565 > #else
566 >   kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
567 > #endif                  
568   }
569  
570  
571  
572 < Photon* find1Photon (PhotonMap *pmap, const RAY* ray)
572 > void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon)
573   {
574 <   float fpos [3], norm [3];
575 <   Photon* photon = NULL;
469 <
470 <   VCOPY(fpos, ray -> rop);
471 <   VCOPY(norm, ray -> ron);
472 <   pmap -> maxDist = thescene.cusize;
473 <   nearest1Neighbour(pmap, fpos, norm, &photon, 1);
474 <  
475 <   return photon;
476 < }
477 <
478 <
479 <
480 < static unsigned long medianPartition (const Photon* heap,
481 <                                      unsigned long* heapIdx,
482 <                                      unsigned long* heapXdi,
483 <                                      unsigned long left,
484 <                                      unsigned long right, unsigned dim)
485 < /* Returns index to median in heap from indices left to right
486 <   (inclusive) in dimension dim. The heap is partitioned relative to
487 <   median using a quicksort algorithm. The heap indices in heapIdx are
488 <   sorted rather than the heap itself. */
489 < {
490 <   register const float* p;
491 <   const unsigned long n = right - left + 1;
492 <   register unsigned long l, r, lg2, n2, m;
493 <   register unsigned d;
494 <  
495 <   /* Round down n to nearest power of 2 */
496 <   for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2);
497 <   n2 = 1 << lg2;
498 <  
499 <   /* Determine median position; this takes into account the fact that
500 <      only the last level in the heap can be partially empty, and that
501 <      it fills from left to right */
502 <   m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1));
503 <  
504 <   while (right > left) {
505 <      /* Pivot node */
506 <      p = heap [heapIdx [right]].pos;
507 <      l = left;
508 <      r = right - 1;
574 > #ifdef PMAP_OOC
575 >   if (OOC_GetPhoton(pmap, idx, photon))
576        
577 <      /* l & r converge, swapping elements out of order with respect to
578 <         pivot node. Identical keys are resolved by cycling through
579 <         dim. The convergence point is then the pivot's position. */
580 <      do {
514 <         while (l <= r) {
515 <            d = dim;
516 <            
517 <            while (heap [heapIdx [l]].pos [d] == p [d]) {
518 <               d = (d + 1) % 3;
519 <              
520 <               if (d == dim) {
521 <                  /* Ignore dupes? */
522 <                  error(WARNING, "duplicate keys in photon heap");
523 <                  l++;
524 <                  break;
525 <               }
526 <            }
527 <            
528 <            if (heap [heapIdx [l]].pos [d] < p [d])
529 <               l++;
530 <            else break;
531 <         }
532 <        
533 <         while (r > l) {
534 <            d = dim;
535 <            
536 <            while (heap [heapIdx [r]].pos [d] == p [d]) {
537 <               d = (d + 1) % 3;
538 <              
539 <               if (d == dim) {
540 <                  /* Ignore dupes? */
541 <                  error(WARNING, "duplicate keys in photon heap");
542 <                  r--;
543 <                  break;
544 <               }
545 <            }
546 <            
547 <            if (heap [heapIdx [r]].pos [d] > p [d])
548 <               r--;
549 <            else break;
550 <         }
551 <        
552 <         /* Swap indices (not the nodes they point to) */
553 <         n2 = heapIdx [l];
554 <         heapIdx [l] = heapIdx [r];
555 <         heapIdx [r] = n2;
556 <         /* Update reverse indices */
557 <         heapXdi [heapIdx [l]] = l;
558 <         heapXdi [n2] = r;
559 <      } while (l < r);
560 <      
561 <      /* Swap indices of convergence and pivot nodes */
562 <      heapIdx [r] = heapIdx [l];
563 <      heapIdx [l] = heapIdx [right];
564 <      heapIdx [right] = n2;
565 <      /* Update reverse indices */
566 <      heapXdi [heapIdx [r]] = r;
567 <      heapXdi [heapIdx [l]] = l;
568 <      heapXdi [n2] = right;
569 <      if (l >= m) right = l - 1;
570 <      if (l <= m) left = l + 1;
571 <   }
572 <  
573 <   /* Once left & right have converged at m, we have found the median */
574 <   return m;
577 > #else
578 >   if (kdT_GetPhoton(pmap, idx, photon))
579 > #endif
580 >      error(INTERNAL, "failed photon lookup");
581   }
582  
583  
584  
585 < void buildHeap (Photon* heap, unsigned long* heapIdx,
580 <                unsigned long* heapXdi, const float min [3],
581 <                const float max [3], unsigned long left,
582 <                unsigned long right, unsigned long root)
583 < /* Recursive part of balancePhotons(..). Builds heap from subarray
584 <   defined by indices left and right. min and max are the minimum resp.
585 <   maximum photon positions in the array. root is the index of the
586 <   current subtree's root, which corresponds to the median's 1-based
587 <   index in the heap. heapIdx are the balanced heap indices. The heap
588 <   is accessed indirectly through these. heapXdi are the reverse indices
589 <   from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */
585 > Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
586   {
587 <   float maxLeft [3], minRight [3];
588 <   Photon rootNode;
589 <   unsigned d;
590 <  
591 <   /* Choose median for dimension with largest spread and partition
596 <      accordingly */
597 <   const float d0 = max [0] - min [0],
598 <               d1 = max [1] - min [1],
599 <               d2 = max [2] - min [2];
600 <   const unsigned char dim = d0 > d1 ? d0 > d2 ? 0 : 2
601 <                                     : d1 > d2 ? 1 : 2;
602 <   const unsigned long median =
603 <      left == right ? left
604 <                    : medianPartition(heap, heapIdx, heapXdi, left, right, dim);
605 <  
606 <   /* Place median at root of current subtree. This consists of swapping
607 <      the median and the root nodes and updating the heap indices */
608 <   memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon));
609 <   memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon));
610 <   setPhotonDiscr(rootNode, dim);
611 <   memcpy(heap + root - 1, &rootNode, sizeof(Photon));
612 <   heapIdx [heapXdi [root - 1]] = heapIdx [median];
613 <   heapXdi [heapIdx [median]] = heapXdi [root - 1];
614 <   heapIdx [median] = root - 1;
615 <   heapXdi [root - 1] = median;
616 <  
617 <   /* Update bounds for left and right subtrees and recurse on them */
618 <   for (d = 0; d <= 2; d++)
619 <      if (d == dim)
620 <         maxLeft [d] = minRight [d] = rootNode.pos [d];
621 <      else {
622 <         maxLeft [d] = max [d];
623 <         minRight [d] = min [d];
624 <      }
625 <      
626 <   if (left < median)
627 <      buildHeap(heap, heapIdx, heapXdi, min, maxLeft,
628 <                left, median - 1, root << 1);
629 <                
630 <   if (right > median)
631 <      buildHeap(heap, heapIdx, heapXdi, minRight, max,
632 <                median + 1, right, (root << 1) + 1);
587 > #ifdef PMAP_OOC
588 >   return OOC_GetNearestPhoton(squeue, idx);
589 > #else
590 >   return kdT_GetNearestPhoton(squeue, idx);
591 > #endif
592   }
593  
594  
595  
596 < void balancePhotons (PhotonMap* pmap, double *photonFlux)
596 > PhotonIdx firstPhoton (const PhotonMap *pmap)
597   {
598 <   Photon *heap = pmap -> heap;
599 <   unsigned long i;
600 <   unsigned long *heapIdx;        /* Photon index array */
601 <   unsigned long *heapXdi;        /* Reverse index to heapIdx */
602 <   unsigned j;
644 <   COLOR flux;
645 <   /* Need doubles here to reduce errors from increment */
646 <   double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
647 <   FVECT d;
648 <  
649 <   if (pmap -> heapEnd) {
650 <      pmap -> heapSize = pmap -> heapEnd;
651 <      heapIdx = (unsigned long*)malloc(pmap -> heapSize *
652 <                                       sizeof(unsigned long));
653 <      heapXdi = (unsigned long*)malloc(pmap -> heapSize *
654 <                                       sizeof(unsigned long));
655 <      if (!heapIdx || !heapXdi)
656 <         error(USER, "can't allocate heap index");
657 <        
658 <      for (i = 0; i < pmap -> heapSize; i++) {
659 <         /* Initialize index arrays */
660 <         heapXdi [i] = heapIdx [i] = i;
661 <         getPhotonFlux(heap + i, flux);
662 <
663 <         /* Scale photon's flux (hitherto normalised to 1 over RGB); in case
664 <          * of a contrib photon map, this is done per light source, and
665 <          * photonFlux is assumed to be an array */
666 <         if (photonFlux) {
667 <            scalecolor(flux, photonFlux [isContribPmap(pmap) ?
668 <                                         photonSrcIdx(pmap, heap + i) : 0]);
669 <            setPhotonFlux(heap + i, flux);
670 <         }
671 <
672 <         /* Need a double here */
673 <         addcolor(avgFlux, flux);
674 <        
675 <         /* Add photon position to centre of gravity */
676 <         for (j = 0; j < 3; j++)
677 <            CoG [j] += heap [i].pos [j];
678 <      }
679 <      
680 <      /* Average photon positions to get centre of gravity */
681 <      for (j = 0; j < 3; j++)
682 <         pmap -> CoG [j] = CoG [j] /= pmap -> heapSize;
683 <      
684 <      /* Compute average photon distance to CoG */
685 <      for (i = 0; i < pmap -> heapSize; i++) {
686 <         VSUB(d, heap [i].pos, CoG);
687 <         CoGdist += DOT(d, d);
688 <      }
689 <      
690 <      pmap -> CoGdist = CoGdist /= pmap -> heapSize;
691 <      
692 <      /* Average photon flux based on RGBE representation */
693 <      scalecolor(avgFlux, 1.0 / pmap -> heapSize);
694 <      copycolor(pmap -> photonFlux, avgFlux);
695 <      
696 <      /* Build kd-tree */
697 <      buildHeap(pmap -> heap, heapIdx, heapXdi, pmap -> minPos,
698 <                pmap -> maxPos, 0, pmap -> heapSize - 1, 1);
699 <                
700 <      free(heapIdx);
701 <      free(heapXdi);
702 <   }
598 > #ifdef PMAP_OOC
599 >   return OOC_FirstPhoton(pmap);
600 > #else
601 >   return kdT_FirstPhoton(pmap);
602 > #endif
603   }
604  
605  
606  
607   void deletePhotons (PhotonMap* pmap)
608   {
609 <   free(pmap -> heap);
610 <   free(pmap -> squeue);
609 > #ifdef PMAP_OOC
610 >   OOC_Delete(&pmap -> store);
611 > #else
612 >   kdT_Delete(&pmap -> store);
613 > #endif
614 >
615 >   free(pmap -> squeue.node);
616     free(pmap -> biasCompHist);
617    
618 <   pmap -> heapSize = 0;
619 <   pmap -> minGather = pmap -> maxGather =
715 <      pmap -> squeueSize = pmap -> squeueEnd = 0;
618 >   pmap -> numPhotons = pmap -> minGather = pmap -> maxGather =
619 >      pmap -> squeue.len = pmap -> squeue.tail = 0;
620   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines