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

Diff Legend

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