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.17 by greg, Thu Sep 29 20:51:07 2016 UTC vs.
Revision 2.23 by rschregle, Wed Apr 14 11:26:25 2021 UTC

# Line 13 | Line 13 | static const char RCSid[] = "$Id$";
13  
14     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
15     (c) Fraunhofer Institute for Solar Energy Systems,
16 +       supported by the German Research Foundation
17 +       (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS))
18     (c) Lucerne University of Applied Sciences and Arts,
19 <       supported by the Swiss National Science Foundation (SNSF, #147053)
19 >       supported by the Swiss National Science Foundation
20 >       (SNSF #147053, "Daylight Redirecting Components")
21     ==========================================================================
22    
23     $Id$
# Line 22 | Line 25 | static const char RCSid[] = "$Id$";
25  
26  
27  
28 < #include "pmap.h"
28 > #include "pmapdata.h"
29   #include "pmaprand.h"
30   #include "pmapmat.h"
31   #include "otypes.h"
# Line 45 | Line 48 | PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
48     #include "pmapkdt.c"
49   #endif
50  
51 + /* Ambient include/exclude set (from ambient.c) */
52 + #ifndef  MAXASET
53 +   #define MAXASET  4095
54 + #endif
55 + extern OBJECT ambset [MAXASET+1];
56  
57 + /* Callback to print photon attributes acc. to user defined format */
58 + int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm);
59  
60 +
61 +
62   void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
63   /* Init photon map 'n' stuff... */
64   {
# Line 70 | Line 82 | void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
82     pmap -> randState [0] = 10243;
83     pmap -> randState [1] = 39829;
84     pmap -> randState [2] = 9433;
73   /* pmapSeed(25999, pmap -> randState); */
85     pmapSeed(randSeed, pmap -> randState);
86    
87     /* Set up type-specific photon lookup callback */
# Line 103 | Line 114 | void initPhotonHeap (PhotonMap *pmap)
114        
115     if (!pmap -> heap) {
116        /* Open heap file */
117 <      if (!(pmap -> heap = tmpfile()))
117 >      mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME));
118 >      if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b")))
119           error(SYSTEM, "failed opening heap file in initPhotonHeap");
120 +
121   #ifdef F_SETFL  /* XXX is there an alternate needed for Windows? */
122        fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
123        fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND);
124 < #endif
112 < /*      ftruncate(fileno(pmap -> heap), 0); */
124 > #endif/*      ftruncate(fileno(pmap -> heap), 0); */
125     }
126   }
127  
# Line 123 | Line 135 | void flushPhotonHeap (PhotonMap *pmap)
135     if (!pmap)
136        error(INTERNAL, "undefined photon map in flushPhotonHeap");
137  
138 <   if (!pmap -> heap || !pmap -> heapBuf)
139 <      error(INTERNAL, "undefined heap in flushPhotonHeap");
138 >   if (!pmap -> heap || !pmap -> heapBuf) {
139 >      /* Silently ignore undefined heap
140 >      error(INTERNAL, "undefined heap in flushPhotonHeap"); */
141 >      return;
142 >   }
143  
144     /* Atomically seek and write block to heap */
145     /* !!! Unbuffered I/O via pwrite() avoids potential race conditions
# Line 142 | Line 157 | void flushPhotonHeap (PhotonMap *pmap)
157     if (write(fd, pmap -> heapBuf, len) != len)
158        error(SYSTEM, "failed append to heap file in flushPhotonHeap");
159  
160 < #if !defined(_WIN32) && !defined(_WIN64)
160 > #if NIX
161     if (fsync(fd))
162        error(SYSTEM, "failed fsync in flushPhotonHeap");
163   #endif
164 +
165     pmap -> heapBufLen = 0;
166   }
167  
168  
169  
170 < #ifdef DEBUG_OOC
170 > #ifdef DEBUG_PMAP
171   static int checkPhotonHeap (FILE *file)
172   /* Check heap for nonsensical or duplicate photons */
173   {
# Line 204 | Line 220 | int newPhoton (PhotonMap* pmap, const RAY* ray)
220     if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
221        return -1;
222  
223 < #ifdef PMAP_ROI
224 <   /* Store photon if within region of interest -- for Ze Eckspertz only! */
225 <   if (ray -> rop [0] >= pmapROI [0] && ray -> rop [0] <= pmapROI [1] &&
226 <       ray -> rop [1] >= pmapROI [2] && ray -> rop [1] <= pmapROI [3] &&
211 <       ray -> rop [2] >= pmapROI [4] && ray -> rop [2] <= pmapROI [5])
212 < #endif
213 <   {      
214 <      /* Adjust flux according to distribution ratio and ray weight */
215 <      copycolor(photonFlux, ray -> rcol);  
216 <      scalecolor(photonFlux,
217 <                 ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
218 <                                                        : 1));
219 <      setPhotonFlux(&photon, photonFlux);
220 <              
221 <      /* Set photon position and flags */
222 <      VCOPY(photon.pos, ray -> rop);
223 <      photon.flags = 0;
224 <      photon.caustic = PMAP_CAUSTICRAY(ray);
223 >   /* Ignore photon if modifier in/outside exclude/include set */
224 >   if (ambincl != -1 && ray -> ro &&
225 >       ambincl != inset(ambset, ray -> ro -> omod))
226 >      return -1;
227  
228 <      /* Set contrib photon's primary ray and subprocess index (the latter
229 <       * to linearise the primary ray indices after photon distribution is
230 <       * complete). Also set primary ray's source index, thereby marking it
229 <       * as used. */
230 <      if (isContribPmap(pmap)) {
231 <         photon.primary = pmap -> numPrimary;
232 <         photon.proc = PMAP_GETRAYPROC(ray);
233 <         pmap -> lastPrimary.srcIdx = ray -> rsrc;
234 <      }
235 <      else photon.primary = 0;
228 >   if (pmapNumROI && pmapROI) {
229 >      unsigned inROI = 0;
230 >      FVECT    photonDist;
231        
232 <      /* Set normal */
233 <      for (i = 0; i <= 2; i++)
234 <         photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
235 <                                                       : ray -> ron [i]);
232 >      /* Store photon if within a region of interest (for ze Ecksperts!)
233 >         Note size of spherical ROI is squared. */
234 >      for (i = 0; !inROI && i < pmapNumROI; i++) {
235 >         VSUB(photonDist, ray -> rop, pmapROI [i].pos);
236 >        
237 >         inROI = (
238 >            PMAP_ROI_ISSPHERE(pmapROI + i)
239 >            ?  DOT(photonDist, photonDist) <= pmapROI [i].siz [0]
240 >            :  fabs(photonDist [0]) <= pmapROI [i].siz [0] &&
241 >               fabs(photonDist [1]) <= pmapROI [i].siz [1] &&
242 >               fabs(photonDist [2]) <= pmapROI [i].siz [2]
243 >         );
244 >      }
245 >      if (!inROI)
246 >         return -1;
247 >   }
248 >  
249 >   /* Adjust flux according to distribution ratio and ray weight */
250 >   copycolor(photonFlux, ray -> rcol);        
251 > #if 0
252 >   /* Factored out ray -> rweight as deprecated (?) for pmap, and infact
253 >      erroneously attenuates volume photon flux based on extinction,
254 >      which is already factored in by photonParticipate() */
255 >   scalecolor(photonFlux,
256 >              ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
257 >                                                     : 1));
258 > #else
259 >   scalecolor(photonFlux,
260 >              1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1));
261 > #endif
262 >   setPhotonFlux(&photon, photonFlux);
263  
264 <      if (!pmap -> heapBuf) {
265 <         /* Lazily allocate heap buffa */
266 < #if 1        
267 <         /* Randomise buffa size to temporally decorellate buffa flushes */        
268 <         srandom(randSeed + getpid());
269 <         pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
264 >   /* Set photon position and flags */
265 >   VCOPY(photon.pos, ray -> rop);
266 >   photon.flags = 0;
267 >   photon.caustic = PMAP_CAUSTICRAY(ray);
268 >
269 >   /* Set contrib photon's primary ray and subprocess index (the latter
270 >    * to linearise the primary ray indices after photon distribution is
271 >    * complete). Also set primary ray's source index, thereby marking it
272 >    * as used. */
273 >   if (isContribPmap(pmap)) {
274 >      photon.primary = pmap -> numPrimary;
275 >      photon.proc = PMAP_GETRAYPROC(ray);
276 >      pmap -> lastPrimary.srcIdx = ray -> rsrc;
277 >   }
278 >   else photon.primary = 0;
279 >  
280 >   /* Set normal */
281 >   for (i = 0; i <= 2; i++)
282 >      photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
283 >                                                    : ray -> ron [i]);
284 >
285 >   if (!pmap -> heapBuf) {
286 >      /* Lazily allocate heap buffa */
287 > #if NIX
288 >      /* Randomise buffa size to temporally decorellate flushes in
289 >       * multiprocessing mode */
290 >      srandom(randSeed + getpid());
291 >      pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
292   #else
293 <         /* Randomisation disabled for reproducability during debugging */        
294 <         pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
293 >      /* Randomisation disabled for single processes on WIN; also useful
294 >       * for reproducability during debugging */        
295 >      pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
296   #endif        
297 <         if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
298 <            error(SYSTEM, "failed heap buffer allocation in newPhoton");
299 <         pmap -> heapBufLen = 0;      
300 <      }
297 >      if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
298 >         error(SYSTEM, "failed heap buffer allocation in newPhoton");
299 >      pmap -> heapBufLen = 0;      
300 >   }
301  
302 <      /* Photon initialised; now append to heap buffa */
303 <      memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
304 <                  
305 <      if (++pmap -> heapBufLen >= pmap -> heapBufSize)
306 <         /* Heap buffa full, flush to heap file */
307 <         flushPhotonHeap(pmap);
302 >   /* Photon initialised; now append to heap buffa */
303 >   memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
304 >              
305 >   if (++pmap -> heapBufLen >= pmap -> heapBufSize)
306 >      /* Heap buffa full, flush to heap file */
307 >      flushPhotonHeap(pmap);
308  
309 <      pmap -> numPhotons++;
310 <   }
309 >   pmap -> numPhotons++;
310 >        
311 >   /* Print photon attributes */
312 >   if (printPhoton)
313 >      /* Non-const kludge */
314 >      printPhoton((RAY*)ray, &photon, pmap);
315              
316     return 0;
317   }
# Line 276 | Line 325 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
325     unsigned       i;
326     Photon         *p;
327     COLOR          flux;
328 +   char           nuHeapFname [sizeof(PMAP_TMPFNAME)];
329     FILE           *nuHeap;
330     /* Need double here to reduce summation errors */
331     double         avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
# Line 285 | Line 335 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
335        error(INTERNAL, "undefined photon map in buildPhotonMap");
336        
337     /* Get number of photons from heapfile size */
338 <   fseek(pmap -> heap, 0, SEEK_END);      
338 >   if (fseek(pmap -> heap, 0, SEEK_END) < 0)
339 >      error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap");
340     pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);
341    
342     if (!pmap -> numPhotons)
# Line 301 | Line 352 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
352     sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons);
353     eputs(errmsg);
354   #endif
355 <    
355 >
356     /* Allocate heap buffa */
357     if (!pmap -> heapBuf) {
358        pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
# Line 313 | Line 364 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
364  
365     /* We REALLY don't need yet another @%&*! heap just to hold the scaled
366      * photons, but can't think of a quicker fix... */
367 <   if (!(nuHeap = tmpfile()))
367 >   mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME));
368 >   if (!(nuHeap = fopen(nuHeapFname, "w+b")))
369        error(SYSTEM, "failed to open postprocessed photon heap in "
370              "buildPhotonMap");
371              
# Line 323 | Line 375 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
375     eputs("Postprocessing photons...\n");
376   #endif
377    
378 <   while (!feof(pmap -> heap)) {
378 >   while (!feof(pmap -> heap)) {  
379 > #ifdef DEBUG_PMAP
380 >      printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap));
381 > #endif      
382        pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
383 <                                 PMAP_HEAPBUFSIZE, pmap -> heap);
384 <      
385 <      if (pmap -> heapBufLen) {
386 <         for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
332 <            /* Update min and max pos and set photon flux */
333 <            for (i = 0; i <= 2; i++) {
334 <               if (p -> pos [i] < pmap -> minPos [i])
335 <                  pmap -> minPos [i] = p -> pos [i];
336 <               else if (p -> pos [i] > pmap -> maxPos [i])
337 <                  pmap -> maxPos [i] = p -> pos [i];  
383 >                                 pmap -> heapBufSize, pmap -> heap);
384 > #ifdef DEBUG_PMAP                                
385 >      printf("Got %lu\n", pmap -> heapBufLen);
386 > #endif      
387  
388 <               /* Update centre of gravity with photon position */                
389 <               CoG [i] += p -> pos [i];                  
341 <            }  
342 <            
343 <            if (primaryOfs)
344 <               /* Linearise photon primary index from subprocess index using the
345 <                * per-subprocess offsets in primaryOfs */
346 <               p -> primary += primaryOfs [p -> proc];
347 <            
348 <            /* Scale photon's flux (hitherto normalised to 1 over RGB); in
349 <             * case of a contrib photon map, this is done per light source,
350 <             * and photonFlux is assumed to be an array */
351 <            getPhotonFlux(p, flux);            
388 >      if (ferror(pmap -> heap))
389 >         error(SYSTEM, "failed to read photon heap in buildPhotonMap");
390  
391 <            if (photonFlux) {
392 <               scalecolor(flux, photonFlux [isContribPmap(pmap) ?
393 <                                               photonSrcIdx(pmap, p) : 0]);
394 <               setPhotonFlux(p, flux);
395 <            }
391 >      for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
392 >         /* Update min and max pos and set photon flux */
393 >         for (i = 0; i <= 2; i++) {
394 >            if (p -> pos [i] < pmap -> minPos [i])
395 >               pmap -> minPos [i] = p -> pos [i];
396 >            else if (p -> pos [i] > pmap -> maxPos [i])
397 >               pmap -> maxPos [i] = p -> pos [i];  
398  
399 <            /* Update average photon flux; need a double here */
400 <            addcolor(avgFlux, flux);
399 >            /* Update centre of gravity with photon position */                
400 >            CoG [i] += p -> pos [i];                  
401 >         }  
402 >        
403 >         if (primaryOfs)
404 >            /* Linearise photon primary index from subprocess index using the
405 >             * per-subprocess offsets in primaryOfs */
406 >            p -> primary += primaryOfs [p -> proc];
407 >        
408 >         /* Scale photon's flux (hitherto normalised to 1 over RGB); in
409 >          * case of a contrib photon map, this is done per light source,
410 >          * and photonFlux is assumed to be an array */
411 >         getPhotonFlux(p, flux);
412 >
413 >         if (photonFlux) {
414 >            scalecolor(flux, photonFlux [isContribPmap(pmap) ?
415 >                                            photonSrcIdx(pmap, p) : 0]);
416 >            setPhotonFlux(p, flux);
417           }
418 +
419 +         /* Update average photon flux; need a double here */
420 +         addcolor(avgFlux, flux);
421 +      }
422          
423 <         /* Write modified photons to new heap */
424 <         fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
423 >      /* Write modified photons to new heap */
424 >      fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
425                  
426 <         if (ferror(nuHeap))
427 <            error(SYSTEM, "failed postprocessing photon flux in "
428 <                  "buildPhotonMap");
369 <      }
426 >      if (ferror(nuHeap))
427 >         error(SYSTEM, "failed postprocessing photon flux in "
428 >               "buildPhotonMap");
429        
430        nCheck += pmap -> heapBufLen;
431     }
# Line 389 | Line 448 | void buildPhotonMap (PhotonMap *pmap, double *photonFl
448     /* Compute average photon distance to centre of gravity */
449     while (!feof(pmap -> heap)) {
450        pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
451 <                                 PMAP_HEAPBUFSIZE, pmap -> heap);
451 >                                 pmap -> heapBufSize, pmap -> heap);
452        
453 <      if (pmap -> heapBufLen)
454 <         for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
455 <            VSUB(d, p -> pos, CoG);
456 <            CoGdist += DOT(d, d);
398 <         }
453 >      for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
454 >         VSUB(d, p -> pos, CoG);
455 >         CoGdist += DOT(d, d);
456 >      }
457     }  
458  
459     pmap -> CoGdist = CoGdist /= pmap -> numPhotons;
460  
461 <   /* Swap heaps */
461 >   /* Swap heaps, discarding unscaled photons */
462     fclose(pmap -> heap);
463 +   unlink(pmap -> heapFname);
464     pmap -> heap = nuHeap;
465 +   strcpy(pmap -> heapFname, nuHeapFname);
466    
467   #ifdef PMAP_OOC
468     OOC_BuildPhotonMap(pmap, nproc);
469   #else
410   /* kd-tree not parallelised */
470     kdT_BuildPhotonMap(pmap);
471   #endif
472  
473     /* Trash heap and its buffa */
474     free(pmap -> heapBuf);
475     fclose(pmap -> heap);
476 +   unlink(pmap -> heapFname);
477     pmap -> heap = NULL;
478     pmap -> heapBuf = NULL;
479   }
# Line 466 | Line 526 | void findPhotons (PhotonMap* pmap, const RAY* ray)
526          
527        /* Search position is ray -> rorg for volume photons, since we have no
528           intersection point. Normals are ignored -- these are incident
529 <         directions). */  
529 >         directions). */
530 >      /* NOTE: status returned by XXX_FindPhotons() is currently ignored;
531 >         if no photons are found, an empty queue is returned under the
532 >         assumption all photons are too distant to contribute significant
533 >         flux. */
534        if (isVolumePmap(pmap)) {
535   #ifdef PMAP_OOC
536           OOC_FindPhotons(pmap, ray -> rorg, NULL);
# Line 556 | Line 620 | void findPhotons (PhotonMap* pmap, const RAY* ray)
620  
621  
622  
623 < void find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
623 > Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
624   {
625 <   pmap -> maxDist2 = thescene.cusize;  /* ? */
625 >   /* Init (squared) search radius to avg photon dist to centre of gravity */
626 >   float maxDist2_0 = pmap -> CoGdist;  
627 >   int found = 0;  
628 > #ifdef PMAP_LOOKUP_REDO
629 >   #define REDO 1
630 > #else
631 >   #define REDO 0
632 > #endif
633    
634 +   do {
635 +      pmap -> maxDist2 = maxDist2_0;
636   #ifdef PMAP_OOC
637 <   OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
637 >      found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
638   #else
639 <   kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
640 < #endif                  
639 >      found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
640 > #endif
641 >      if (found < 0) {
642 >         /* Expand search radius to retry */
643 >         maxDist2_0 *= 2;      
644 > #ifdef PMAP_LOOKUP_WARN
645 >         sprintf(errmsg, "failed 1-NN photon lookup"
646 > #ifdef PMAP_LOOKUP_REDO
647 >            ", retrying with search radius %.2f", maxDist2_0
648 > #endif
649 >         );
650 >         error(WARNING, errmsg);
651 > #endif
652 >      }
653 >   } while (REDO && found < 0);
654 >
655 >   /* Return photon buffer containing valid photon, else NULL */
656 >   return found < 0 ? NULL : photon;
657   }
658  
659  
# Line 572 | Line 661 | void find1Photon (PhotonMap *pmap, const RAY* ray, Pho
661   void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon)
662   {
663   #ifdef PMAP_OOC
664 <   if (OOC_GetPhoton(pmap, idx, photon))
576 <      
664 >   if (OOC_GetPhoton(pmap, idx, photon))      
665   #else
666     if (kdT_GetPhoton(pmap, idx, photon))
667   #endif

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines