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.15 by rschregle, Tue May 17 17:39:47 2016 UTC vs.
Revision 2.22 by rschregle, Wed Apr 8 15:14:21 2020 UTC

# Line 22 | Line 22 | static const char RCSid[] = "$Id$";
22  
23  
24  
25 < #include "pmap.h"
25 > #include "pmapdata.h"
26   #include "pmaprand.h"
27   #include "pmapmat.h"
28   #include "otypes.h"
# Line 45 | Line 45 | PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
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 + /* Callback to print photon attributes acc. to user defined format */
55 + int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm);
56  
57 +
58 +
59   void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
60   /* Init photon map 'n' stuff... */
61   {
# Line 70 | Line 79 | void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
79     pmap -> randState [0] = 10243;
80     pmap -> randState [1] = 39829;
81     pmap -> randState [2] = 9433;
73   /* pmapSeed(25999, pmap -> randState); */
82     pmapSeed(randSeed, pmap -> randState);
83    
84     /* Set up type-specific photon lookup callback */
# Line 103 | Line 111 | void initPhotonHeap (PhotonMap *pmap)
111        
112     if (!pmap -> heap) {
113        /* Open heap file */
114 <      if (!(pmap -> heap = tmpfile()))
114 >      mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME));
115 >      if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b")))
116           error(SYSTEM, "failed opening heap file in initPhotonHeap");
117 +
118 + #ifdef F_SETFL  /* XXX is there an alternate needed for Windows? */
119        fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
120        fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND);
121 < /*      ftruncate(fileno(pmap -> heap), 0); */
121 > #endif/*      ftruncate(fileno(pmap -> heap), 0); */
122     }
123   }
124  
# Line 121 | Line 132 | void flushPhotonHeap (PhotonMap *pmap)
132     if (!pmap)
133        error(INTERNAL, "undefined photon map in flushPhotonHeap");
134  
135 <   if (!pmap -> heap || !pmap -> heapBuf)
136 <      error(INTERNAL, "undefined heap in flushPhotonHeap");
135 >   if (!pmap -> heap || !pmap -> heapBuf) {
136 >      /* Silently ignore undefined heap
137 >      error(INTERNAL, "undefined heap in flushPhotonHeap"); */
138 >      return;
139 >   }
140  
141     /* Atomically seek and write block to heap */
142     /* !!! Unbuffered I/O via pwrite() avoids potential race conditions
# Line 139 | Line 153 | void flushPhotonHeap (PhotonMap *pmap)
153     /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */
154     if (write(fd, pmap -> heapBuf, len) != len)
155        error(SYSTEM, "failed append to heap file in flushPhotonHeap");
156 <  
156 >
157 > #if NIX
158     if (fsync(fd))
159        error(SYSTEM, "failed fsync in flushPhotonHeap");
160 <      
160 > #endif
161 >
162     pmap -> heapBufLen = 0;
163   }
164  
165  
166  
167 < #ifdef DEBUG_OOC
167 > #ifdef DEBUG_PMAP
168   static int checkPhotonHeap (FILE *file)
169   /* Check heap for nonsensical or duplicate photons */
170   {
# Line 201 | Line 217 | int newPhoton (PhotonMap* pmap, const RAY* ray)
217     if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
218        return -1;
219  
220 < #ifdef PMAP_ROI
221 <   /* Store photon if within region of interest -- for Ze Eckspertz only! */
222 <   if (ray -> rop [0] >= pmapROI [0] && ray -> rop [0] <= pmapROI [1] &&
223 <       ray -> rop [1] >= pmapROI [2] && ray -> rop [1] <= pmapROI [3] &&
208 <       ray -> rop [2] >= pmapROI [4] && ray -> rop [2] <= pmapROI [5])
209 < #endif
210 <   {      
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);
217 <              
218 <      /* Set photon position and flags */
219 <      VCOPY(photon.pos, ray -> rop);
220 <      photon.flags = 0;
221 <      photon.caustic = PMAP_CAUSTICRAY(ray);
220 >   /* Ignore photon if modifier in/outside exclude/include set */
221 >   if (ambincl != -1 && ray -> ro &&
222 >       ambincl != inset(ambset, ray -> ro -> omod))
223 >      return -1;
224  
225 <      /* Set contrib photon's primary ray and subprocess index (the latter
226 <       * 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 -> numPrimary;
229 <         photon.proc = PMAP_GETRAYPROC(ray);
230 <         pmap -> lastPrimary.srcIdx = ray -> rsrc;
231 <      }
232 <      else photon.primary = 0;
225 >   if (pmapNumROI && pmapROI) {      
226 >      unsigned inROI = 0;
227        
228 <      /* Set normal */
229 <      for (i = 0; i <= 2; i++)
230 <         photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
231 <                                                       : ray -> ron [i]);
228 >      /* Store photon if within a region of interest (for ze Ecksperts!) */
229 >      for (i = 0; !inROI && i < pmapNumROI; i++)
230 >         inROI = (ray -> rop [0] >= pmapROI [i].min [0] &&
231 >                  ray -> rop [0] <= pmapROI [i].max [0] &&
232 >                  ray -> rop [1] >= pmapROI [i].min [1] &&
233 >                  ray -> rop [1] <= pmapROI [i].max [1] &&
234 >                  ray -> rop [2] >= pmapROI [i].min [2] &&
235 >                  ray -> rop [2] <= pmapROI [i].max [2]);
236 >      if (!inROI)
237 >         return -1;
238 >   }
239 >  
240 >   /* Adjust flux according to distribution ratio and ray weight */
241 >   copycolor(photonFlux, ray -> rcol);        
242 > #if 0
243 >   /* Factored out ray -> rweight as deprecated (?) for pmap, and infact
244 >      erroneously attenuates volume photon flux based on extinction,
245 >      which is already factored in by photonParticipate() */
246 >   scalecolor(photonFlux,
247 >              ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
248 >                                                     : 1));
249 > #else
250 >   scalecolor(photonFlux,
251 >              1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1));
252 > #endif
253 >   setPhotonFlux(&photon, photonFlux);
254  
255 <      if (!pmap -> heapBuf) {
256 <         /* Lazily allocate heap buffa */
257 < #if 1        
258 <         /* Randomise buffa size to temporally decorellate buffa flushes */        
259 <         srandom(randSeed + getpid());
260 <         pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
255 >   /* Set photon position and flags */
256 >   VCOPY(photon.pos, ray -> rop);
257 >   photon.flags = 0;
258 >   photon.caustic = PMAP_CAUSTICRAY(ray);
259 >
260 >   /* Set contrib photon's primary ray and subprocess index (the latter
261 >    * to linearise the primary ray indices after photon distribution is
262 >    * complete). Also set primary ray's source index, thereby marking it
263 >    * as used. */
264 >   if (isContribPmap(pmap)) {
265 >      photon.primary = pmap -> numPrimary;
266 >      photon.proc = PMAP_GETRAYPROC(ray);
267 >      pmap -> lastPrimary.srcIdx = ray -> rsrc;
268 >   }
269 >   else photon.primary = 0;
270 >  
271 >   /* Set normal */
272 >   for (i = 0; i <= 2; i++)
273 >      photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
274 >                                                    : ray -> ron [i]);
275 >
276 >   if (!pmap -> heapBuf) {
277 >      /* Lazily allocate heap buffa */
278 > #if NIX
279 >      /* Randomise buffa size to temporally decorellate flushes in
280 >       * multiprocessing mode */
281 >      srandom(randSeed + getpid());
282 >      pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
283   #else
284 <         /* Randomisation disabled for reproducability during debugging */        
285 <         pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
284 >      /* Randomisation disabled for single processes on WIN; also useful
285 >       * for reproducability during debugging */        
286 >      pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
287   #endif        
288 <         if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
289 <            error(SYSTEM, "failed heap buffer allocation in newPhoton");
290 <         pmap -> heapBufLen = 0;      
291 <      }
288 >      if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
289 >         error(SYSTEM, "failed heap buffer allocation in newPhoton");
290 >      pmap -> heapBufLen = 0;      
291 >   }
292  
293 <      /* Photon initialised; now append to heap buffa */
294 <      memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
295 <                  
296 <      if (++pmap -> heapBufLen >= pmap -> heapBufSize)
297 <         /* Heap buffa full, flush to heap file */
298 <         flushPhotonHeap(pmap);
293 >   /* Photon initialised; now append to heap buffa */
294 >   memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
295 >              
296 >   if (++pmap -> heapBufLen >= pmap -> heapBufSize)
297 >      /* Heap buffa full, flush to heap file */
298 >      flushPhotonHeap(pmap);
299  
300 <      pmap -> numPhotons++;
301 <   }
300 >   pmap -> numPhotons++;
301 >        
302 >   /* Print photon attributes */
303 >   if (printPhoton)
304 >      /* Non-const kludge */
305 >      printPhoton((RAY*)ray, &photon, pmap);
306              
307     return 0;
308   }
# Line 273 | Line 316 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
316     unsigned       i;
317     Photon         *p;
318     COLOR          flux;
319 +   char           nuHeapFname [sizeof(PMAP_TMPFNAME)];
320     FILE           *nuHeap;
321     /* Need double here to reduce summation errors */
322     double         avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
# Line 282 | Line 326 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
326        error(INTERNAL, "undefined photon map in buildPhotonMap");
327        
328     /* Get number of photons from heapfile size */
329 <   fseek(pmap -> heap, 0, SEEK_END);      
329 >   if (fseek(pmap -> heap, 0, SEEK_END) < 0)
330 >      error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap");
331     pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);
332    
333     if (!pmap -> numPhotons)
# Line 298 | Line 343 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
343     sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons);
344     eputs(errmsg);
345   #endif
346 <    
346 >
347     /* Allocate heap buffa */
348     if (!pmap -> heapBuf) {
349        pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
# Line 310 | Line 355 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
355  
356     /* We REALLY don't need yet another @%&*! heap just to hold the scaled
357      * photons, but can't think of a quicker fix... */
358 <   if (!(nuHeap = tmpfile()))
358 >   mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME));
359 >   if (!(nuHeap = fopen(nuHeapFname, "w+b")))
360        error(SYSTEM, "failed to open postprocessed photon heap in "
361              "buildPhotonMap");
362              
# Line 320 | Line 366 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
366     eputs("Postprocessing photons...\n");
367   #endif
368    
369 <   while (!feof(pmap -> heap)) {
369 >   while (!feof(pmap -> heap)) {  
370 > #ifdef DEBUG_PMAP
371 >      printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap));
372 > #endif      
373        pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
374 <                                 PMAP_HEAPBUFSIZE, pmap -> heap);
375 <      
376 <      if (pmap -> heapBufLen) {
377 <         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];  
374 >                                 pmap -> heapBufSize, pmap -> heap);
375 > #ifdef DEBUG_PMAP                                
376 >      printf("Got %lu\n", pmap -> heapBufLen);
377 > #endif      
378  
379 <               /* Update centre of gravity with photon position */                
380 <               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);            
379 >      if (ferror(pmap -> heap))
380 >         error(SYSTEM, "failed to read photon heap in buildPhotonMap");
381  
382 <            if (photonFlux) {
383 <               scalecolor(flux, photonFlux [isContribPmap(pmap) ?
384 <                                               photonSrcIdx(pmap, p) : 0]);
385 <               setPhotonFlux(p, flux);
386 <            }
382 >      for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
383 >         /* Update min and max pos and set photon flux */
384 >         for (i = 0; i <= 2; i++) {
385 >            if (p -> pos [i] < pmap -> minPos [i])
386 >               pmap -> minPos [i] = p -> pos [i];
387 >            else if (p -> pos [i] > pmap -> maxPos [i])
388 >               pmap -> maxPos [i] = p -> pos [i];  
389  
390 <            /* Update average photon flux; need a double here */
391 <            addcolor(avgFlux, flux);
390 >            /* Update centre of gravity with photon position */                
391 >            CoG [i] += p -> pos [i];                  
392 >         }  
393 >        
394 >         if (primaryOfs)
395 >            /* Linearise photon primary index from subprocess index using the
396 >             * per-subprocess offsets in primaryOfs */
397 >            p -> primary += primaryOfs [p -> proc];
398 >        
399 >         /* Scale photon's flux (hitherto normalised to 1 over RGB); in
400 >          * case of a contrib photon map, this is done per light source,
401 >          * and photonFlux is assumed to be an array */
402 >         getPhotonFlux(p, flux);
403 >
404 >         if (photonFlux) {
405 >            scalecolor(flux, photonFlux [isContribPmap(pmap) ?
406 >                                            photonSrcIdx(pmap, p) : 0]);
407 >            setPhotonFlux(p, flux);
408           }
409 +
410 +         /* Update average photon flux; need a double here */
411 +         addcolor(avgFlux, flux);
412 +      }
413          
414 <         /* Write modified photons to new heap */
415 <         fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
414 >      /* Write modified photons to new heap */
415 >      fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
416                  
417 <         if (ferror(nuHeap))
418 <            error(SYSTEM, "failed postprocessing photon flux in "
419 <                  "buildPhotonMap");
366 <      }
417 >      if (ferror(nuHeap))
418 >         error(SYSTEM, "failed postprocessing photon flux in "
419 >               "buildPhotonMap");
420        
421        nCheck += pmap -> heapBufLen;
422     }
# Line 386 | Line 439 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
439     /* Compute average photon distance to centre of gravity */
440     while (!feof(pmap -> heap)) {
441        pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
442 <                                 PMAP_HEAPBUFSIZE, pmap -> heap);
442 >                                 pmap -> heapBufSize, pmap -> heap);
443        
444 <      if (pmap -> heapBufLen)
445 <         for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
446 <            VSUB(d, p -> pos, CoG);
447 <            CoGdist += DOT(d, d);
395 <         }
444 >      for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
445 >         VSUB(d, p -> pos, CoG);
446 >         CoGdist += DOT(d, d);
447 >      }
448     }  
449  
450     pmap -> CoGdist = CoGdist /= pmap -> numPhotons;
451  
452 <   /* Swap heaps */
452 >   /* Swap heaps, discarding unscaled photons */
453     fclose(pmap -> heap);
454 +   unlink(pmap -> heapFname);
455     pmap -> heap = nuHeap;
456 +   strcpy(pmap -> heapFname, nuHeapFname);
457    
458   #ifdef PMAP_OOC
459     OOC_BuildPhotonMap(pmap, nproc);
460   #else
407   /* kd-tree not parallelised */
461     kdT_BuildPhotonMap(pmap);
462   #endif
463  
464     /* Trash heap and its buffa */
465     free(pmap -> heapBuf);
466     fclose(pmap -> heap);
467 +   unlink(pmap -> heapFname);
468     pmap -> heap = NULL;
469     pmap -> heapBuf = NULL;
470   }
# Line 463 | Line 517 | void findPhotons (PhotonMap* pmap, const RAY* ray)
517          
518        /* Search position is ray -> rorg for volume photons, since we have no
519           intersection point. Normals are ignored -- these are incident
520 <         directions). */  
520 >         directions). */
521 >      /* NOTE: status returned by XXX_FindPhotons() is currently ignored;
522 >         if no photons are found, an empty queue is returned under the
523 >         assumption all photons are too distant to contribute significant
524 >         flux. */
525        if (isVolumePmap(pmap)) {
526   #ifdef PMAP_OOC
527           OOC_FindPhotons(pmap, ray -> rorg, NULL);
# Line 553 | Line 611 | void findPhotons (PhotonMap* pmap, const RAY* ray)
611  
612  
613  
614 < void find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
614 > Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
615   {
616 <   pmap -> maxDist2 = thescene.cusize;  /* ? */
616 >   /* Init (squared) search radius to avg photon dist to centre of gravity */
617 >   float maxDist2_0 = pmap -> CoGdist;  
618 >   int found = 0;  
619 > #ifdef PMAP_LOOKUP_REDO
620 >   #define REDO 1
621 > #else
622 >   #define REDO 0
623 > #endif
624    
625 +   do {
626 +      pmap -> maxDist2 = maxDist2_0;
627   #ifdef PMAP_OOC
628 <   OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
628 >      found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
629   #else
630 <   kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
631 < #endif                  
630 >      found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
631 > #endif
632 >      if (found < 0) {
633 >         /* Expand search radius to retry */
634 >         maxDist2_0 *= 2;      
635 > #ifdef PMAP_LOOKUP_WARN
636 >         sprintf(errmsg, "failed 1-NN photon lookup"
637 > #ifdef PMAP_LOOKUP_REDO
638 >            ", retrying with search radius %.2f", maxDist2_0
639 > #endif
640 >         );
641 >         error(WARNING, errmsg);
642 > #endif
643 >      }
644 >   } while (REDO && found < 0);
645 >
646 >   /* Return photon buffer containing valid photon, else NULL */
647 >   return found < 0 ? NULL : photon;
648   }
649  
650  
# Line 569 | Line 652 | void find1Photon (PhotonMap *pmap, const RAY* ray, Pho
652   void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon)
653   {
654   #ifdef PMAP_OOC
655 <   if (OOC_GetPhoton(pmap, idx, photon))
573 <      
655 >   if (OOC_GetPhoton(pmap, idx, photon))      
656   #else
657     if (kdT_GetPhoton(pmap, idx, photon))
658   #endif

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines