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

Diff Legend

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