src/rt/pmapdata.c

#ifndef lint
static const char RCSid[] = "$Id: pmapdata.c,v 2.17 2016/09/29 20:51:07 greg Exp $";
#endif

/* 
   =========================================================================
   Photon map types and interface to nearest neighbour lookups in underlying
   point cloud data structure.
   
   The default data structure is an in-core kd-tree (see pmapkdt.{h,c}).
   This can be overriden with the PMAP_OOC compiletime switch, which enables
   an out-of-core octree (see oococt.{h,c}).

   Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
   (c) Fraunhofer Institute for Solar Energy Systems,
   (c) Lucerne University of Applied Sciences and Arts,
       supported by the Swiss National Science Foundation (SNSF, #147053)
   ==========================================================================
   
   $Id: pmapdata.c,v 2.17 2016/09/29 20:51:07 greg Exp $
*/


#include "pmapdata.h"
#include "pmaprand.h"
#include "pmapmat.h"
#include "otypes.h"
#include "source.h"
#include "rcontrib.h"
#include "random.h"


PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
   NULL, NULL, NULL, NULL, NULL, NULL
};


/* Include routines to handle underlying point cloud data structure */
#ifdef PMAP_OOC
   #include "pmapooc.c"
#else
   #include "pmapkdt.c"
#endif


void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
/* Init photon map 'n' stuff... */
{
   if (!pmap) 
      return;
      
   pmap -> numPhotons = 0;
   pmap -> biasCompHist = NULL;
   pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE;
   pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE;
   pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0;
   pmap -> gatherTolerance = gatherTolerance;
   pmap -> minError = pmap -> maxError = pmap -> rmsError = 0;
   pmap -> numDensity = 0;
   pmap -> distribRatio = 1;
   pmap -> type = t;
   pmap -> squeue.node = NULL;
   pmap -> squeue.len = 0;   

   /* Init local RNG state */
   pmap -> randState [0] = 10243;
   pmap -> randState [1] = 39829;
   pmap -> randState [2] = 9433;
   pmapSeed(randSeed, pmap -> randState);
   
   /* Set up type-specific photon lookup callback */
   pmap -> lookup = pmapLookup [t];

   /* Mark primary photon ray as unused */
   pmap -> lastPrimary.srcIdx = -1;
   pmap -> numPrimary = 0;  
   pmap -> primaries = NULL;
   
   /* Init storage */
   pmap -> heap = NULL;
   pmap -> heapBuf = NULL;
   pmap -> heapBufLen = 0;
#ifdef PMAP_OOC
   OOC_Null(&pmap -> store);
#else
   kdT_Null(&pmap -> store);
#endif
}


void initPhotonHeap (PhotonMap *pmap)
{
   int fdFlags;
   
   if (!pmap)
      error(INTERNAL, "undefined photon map in initPhotonHeap");
      
   if (!pmap -> heap) {
      /* Open heap file */
      mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME));
      if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b")))
         error(SYSTEM, "failed opening heap file in initPhotonHeap");

#ifdef F_SETFL  /* XXX is there an alternate needed for Windows? */
      fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
      fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND);
#endif/*      ftruncate(fileno(pmap -> heap), 0); */
   }
}


void flushPhotonHeap (PhotonMap *pmap)
{
   int                  fd;
   const unsigned long  len = pmap -> heapBufLen * sizeof(Photon);
   
   if (!pmap)
      error(INTERNAL, "undefined photon map in flushPhotonHeap");

   if (!pmap -> heap || !pmap -> heapBuf) {
      /* Silently ignore undefined heap 
      error(INTERNAL, "undefined heap in flushPhotonHeap"); */
      return;
   }

   /* Atomically seek and write block to heap */
   /* !!! Unbuffered I/O via pwrite() avoids potential race conditions
    * !!! and buffer corruption which can occur with lseek()/fseek()
    * !!! followed by write()/fwrite(). */   
   fd = fileno(pmap -> heap);
   
#ifdef DEBUG_PMAP
   sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), 
           pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon)); 
   eputs(errmsg);
#endif

   /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */
   if (write(fd, pmap -> heapBuf, len) != len)
      error(SYSTEM, "failed append to heap file in flushPhotonHeap");

#if NIX
   if (fsync(fd))
      error(SYSTEM, "failed fsync in flushPhotonHeap");
#endif

   pmap -> heapBufLen = 0;
}


#ifdef DEBUG_PMAP
static int checkPhotonHeap (FILE *file)
/* Check heap for nonsensical or duplicate photons */
{
   Photon   p, lastp;
   int      i, dup;
   
   rewind(file);
   memset(&lastp, 0, sizeof(lastp));
   
   while (fread(&p, sizeof(p), 1, file)) {
      dup = 1;
      
      for (i = 0; i <= 2; i++) {
         if (p.pos [i] < thescene.cuorg [i] || 
             p.pos [i] > thescene.cuorg [i] + thescene.cusize) { 
             
            sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", 
                    p.pos [0], p.pos [1], p.pos [2]);
            error(WARNING, errmsg);
         }
         
         dup &= p.pos [i] == lastp.pos [i];
      }
      
      if (dup) {
         sprintf(errmsg,
                 "consecutive duplicate photon in heap at [%f, %f, %f]\n", 
                 p.pos [0], p.pos [1], p.pos [2]);
         error(WARNING, errmsg);
      }
   }

   return 0;
}
#endif


int newPhoton (PhotonMap* pmap, const RAY* ray)
{
   unsigned i, inROI = 0;
   Photon photon;
   COLOR photonFlux;
   
   /* Account for distribution ratio */
   if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) 
      return -1;
      
   /* Don't store on sources */
   if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) 
      return -1;

   /* Store photon if within a region of interest (for ze Ecksperts!) */
   if (!pmapNumROI || !pmapROI) 
      inROI = 1;
   else {
      for (i = 0; !inROI && i < pmapNumROI; i++)
         inROI = (ray -> rop [0] >= pmapROI [i].min [0] && 
                  ray -> rop [0] <= pmapROI [i].max [0] &&
                  ray -> rop [1] >= pmapROI [i].min [1] && 
                  ray -> rop [1] <= pmapROI [i].max [1] &&
                  ray -> rop [2] >= pmapROI [i].min [2] && 
                  ray -> rop [2] <= pmapROI [i].max [2]);
   }
    
   if (inROI) {       
      /* Adjust flux according to distribution ratio and ray weight */
      copycolor(photonFlux, ray -> rcol);   
      scalecolor(photonFlux, 
                 ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
                                                        : 1));
      setPhotonFlux(&photon, photonFlux);
               
      /* Set photon position and flags */
      VCOPY(photon.pos, ray -> rop);
      photon.flags = 0;
      photon.caustic = PMAP_CAUSTICRAY(ray);

      /* Set contrib photon's primary ray and subprocess index (the latter
       * to linearise the primary ray indices after photon distribution is
       * complete). Also set primary ray's source index, thereby marking it
       * as used. */
      if (isContribPmap(pmap)) {
         photon.primary = pmap -> numPrimary;
         photon.proc = PMAP_GETRAYPROC(ray);
         pmap -> lastPrimary.srcIdx = ray -> rsrc;
      }
      else photon.primary = 0;
      
      /* Set normal */
      for (i = 0; i <= 2; i++)
         photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i] 
                                                       : ray -> ron [i]);

      if (!pmap -> heapBuf) {
         /* Lazily allocate heap buffa */
#if NIX
         /* Randomise buffa size to temporally decorellate flushes in
          * multiprocessing mode */
         srandom(randSeed + getpid());
         pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
#else
         /* Randomisation disabled for single processes on WIN; also useful
          * for reproducability during debugging */         
         pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
#endif         
         if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
            error(SYSTEM, "failed heap buffer allocation in newPhoton");
         pmap -> heapBufLen = 0;      
      }

      /* Photon initialised; now append to heap buffa */
      memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
                  
      if (++pmap -> heapBufLen >= pmap -> heapBufSize)
         /* Heap buffa full, flush to heap file */
         flushPhotonHeap(pmap);

      pmap -> numPhotons++;
   }
            
   return 0;
}


void buildPhotonMap (PhotonMap *pmap, double *photonFlux, 
                     PhotonPrimaryIdx *primaryOfs, unsigned nproc)
{
   unsigned long  n, nCheck = 0;
   unsigned       i;
   Photon         *p;
   COLOR          flux;
   char           nuHeapFname [sizeof(PMAP_TMPFNAME)];
   FILE           *nuHeap;
   /* Need double here to reduce summation errors */
   double         avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
   FVECT          d;
   
   if (!pmap)
      error(INTERNAL, "undefined photon map in buildPhotonMap");
      
   /* Get number of photons from heapfile size */
   if (fseek(pmap -> heap, 0, SEEK_END) < 0)
      error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap");
   pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);
   
   if (!pmap -> numPhotons)
      error(INTERNAL, "empty photon map in buildPhotonMap");   

   if (!pmap -> heap)
      error(INTERNAL, "no heap in buildPhotonMap");

#ifdef DEBUG_PMAP
   eputs("Checking photon heap consistency...\n");
   checkPhotonHeap(pmap -> heap);
   
   sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons);
   eputs(errmsg);
#endif

   /* Allocate heap buffa */
   if (!pmap -> heapBuf) {
      pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
      pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon));
      if (!pmap -> heapBuf)
         error(SYSTEM, "failed to allocate postprocessed photon heap in" 
               "buildPhotonMap");
   }

   /* We REALLY don't need yet another @%&*! heap just to hold the scaled
    * photons, but can't think of a quicker fix... */
   mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME));
   if (!(nuHeap = fopen(nuHeapFname, "w+b")))
      error(SYSTEM, "failed to open postprocessed photon heap in "
            "buildPhotonMap");
            
   rewind(pmap -> heap);

#ifdef DEBUG_PMAP 
   eputs("Postprocessing photons...\n");
#endif
   
   while (!feof(pmap -> heap)) {   
#ifdef DEBUG_PMAP 
      printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap));
#endif      
      pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), 
                                 pmap -> heapBufSize, pmap -> heap);
#ifdef DEBUG_PMAP                                 
      printf("Got %lu\n", pmap -> heapBufLen);
#endif      

      if (ferror(pmap -> heap))
         error(SYSTEM, "failed to read photon heap in buildPhotonMap");

      for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
         /* Update min and max pos and set photon flux */
         for (i = 0; i <= 2; i++) {
            if (p -> pos [i] < pmap -> minPos [i]) 
               pmap -> minPos [i] = p -> pos [i];
            else if (p -> pos [i] > pmap -> maxPos [i]) 
               pmap -> maxPos [i] = p -> pos [i];   

            /* Update centre of gravity with photon position */                 
            CoG [i] += p -> pos [i];                  
         }  
         
         if (primaryOfs)
            /* Linearise photon primary index from subprocess index using the
             * per-subprocess offsets in primaryOfs */
            p -> primary += primaryOfs [p -> proc];
         
         /* Scale photon's flux (hitherto normalised to 1 over RGB); in
          * case of a contrib photon map, this is done per light source,
          * and photonFlux is assumed to be an array */
         getPhotonFlux(p, flux);            

         if (photonFlux) {
            scalecolor(flux, photonFlux [isContribPmap(pmap) ? 
                                            photonSrcIdx(pmap, p) : 0]);
            setPhotonFlux(p, flux);
         }

         /* Update average photon flux; need a double here */
         addcolor(avgFlux, flux);
      }
         
      /* Write modified photons to new heap */
      fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
                
      if (ferror(nuHeap))
         error(SYSTEM, "failed postprocessing photon flux in "
               "buildPhotonMap");
      
      nCheck += pmap -> heapBufLen;
   }

#ifdef DEBUG_PMAP
   if (nCheck < pmap -> numPhotons)
      error(INTERNAL, "truncated photon heap in buildPhotonMap");
#endif
   
   /* Finalise average photon flux */
   scalecolor(avgFlux, 1.0 / pmap -> numPhotons);
   copycolor(pmap -> photonFlux, avgFlux);

   /* Average photon positions to get centre of gravity */
   for (i = 0; i < 3; i++)
      pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons;
      
   rewind(pmap -> heap);
   
   /* Compute average photon distance to centre of gravity */
   while (!feof(pmap -> heap)) {
      pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), 
                                 pmap -> heapBufSize, pmap -> heap);
      
      for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
         VSUB(d, p -> pos, CoG);
         CoGdist += DOT(d, d);
      }
   }   

   pmap -> CoGdist = CoGdist /= pmap -> numPhotons;

   /* Swap heaps, discarding unscaled photons */
   fclose(pmap -> heap);
   unlink(pmap -> heapFname);
   pmap -> heap = nuHeap;
   strcpy(pmap -> heapFname, nuHeapFname);
   
#ifdef PMAP_OOC
   OOC_BuildPhotonMap(pmap, nproc);
#else
   kdT_BuildPhotonMap(pmap);
#endif

   /* Trash heap and its buffa */
   free(pmap -> heapBuf);
   fclose(pmap -> heap);
   unlink(pmap -> heapFname);
   pmap -> heap = NULL;
   pmap -> heapBuf = NULL;
}


/* Dynamic max photon search radius increase and reduction factors */
#define PMAP_MAXDIST_INC 4
#define PMAP_MAXDIST_DEC 0.9

/* Num successful lookups before reducing in max search radius */
#define PMAP_MAXDIST_CNT 1000

/* Threshold below which we assume increasing max radius won't help */
#define PMAP_SHORT_LOOKUP_THRESH 1

/* Coefficient for adaptive maximum search radius */
#define PMAP_MAXDIST_COEFF 100

void findPhotons (PhotonMap* pmap, const RAY* ray)
{
   int redo = 0;
   
   if (!pmap -> squeue.len) {
      /* Lazy init priority queue */
#ifdef PMAP_OOC
      OOC_InitFindPhotons(pmap);
#else
      kdT_InitFindPhotons(pmap);      
#endif      
      pmap -> minGathered = pmap -> maxGather;
      pmap -> maxGathered = pmap -> minGather;
      pmap -> totalGathered = 0;
      pmap -> numLookups = pmap -> numShortLookups = 0;
      pmap -> shortLookupPct = 0;
      pmap -> minError = FHUGE;
      pmap -> maxError = -FHUGE;
      pmap -> rmsError = 0;
      /* SQUARED max search radius limit is based on avg photon distance to
       * centre of gravity, unless fixed by user (maxDistFix > 0) */
      pmap -> maxDist0 = pmap -> maxDist2Limit = 
         maxDistFix > 0 ? maxDistFix * maxDistFix
                        : PMAP_MAXDIST_COEFF * pmap -> squeue.len * 
                          pmap -> CoGdist / pmap -> numPhotons;
   }

   do {
      pmap -> squeue.tail = 0;
      pmap -> maxDist2 = pmap -> maxDist0;
         
      /* Search position is ray -> rorg for volume photons, since we have no 
         intersection point. Normals are ignored -- these are incident 
         directions). */   
      if (isVolumePmap(pmap)) {
#ifdef PMAP_OOC
         OOC_FindPhotons(pmap, ray -> rorg, NULL);
#else
         kdT_FindPhotons(pmap, ray -> rorg, NULL);
#endif                  
      }
      else {
#ifdef PMAP_OOC
         OOC_FindPhotons(pmap, ray -> rop, ray -> ron);         
#else
         kdT_FindPhotons(pmap, ray -> rop, ray -> ron);
#endif
      }

#ifdef PMAP_LOOKUP_INFO
      fprintf(stderr, "%d/%d %s photons found within radius %.3f " 
              "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, 
              pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2),
              ray -> rop [0], ray -> rop [1], ray -> rop [2], 
              ray -> ro ? ray -> ro -> oname : "<null>");
#endif      
            
      if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) {
         /* Short lookup; too few photons found */
         if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) {
            /* Ignore short lookups which return fewer than
             * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there
             * really are no photons in the vicinity, and increasing the max
             * search radius therefore won't help */
#ifdef PMAP_LOOKUP_WARN
            sprintf(errmsg, 
                    "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
                    pmap -> squeue.tail, pmap -> squeue.len, 
                    pmapName [pmap -> type], 
                    ray -> rop [0], ray -> rop [1], ray -> rop [2], 
                    ray -> ro ? ray -> ro -> oname : "<null>");
            error(WARNING, errmsg);         
#endif             

            /* Bail out after warning if maxDist is fixed */
            if (maxDistFix > 0)
               return;
            
            if (pmap -> maxDist0 < pmap -> maxDist2Limit) {
               /* Increase max search radius if below limit & redo search */
               pmap -> maxDist0 *= PMAP_MAXDIST_INC;
#ifdef PMAP_LOOKUP_REDO
               redo = 1;
#endif               
#ifdef PMAP_LOOKUP_WARN
               sprintf(errmsg, 
                       redo ? "restarting photon lookup with max radius %.1e"
                            : "max photon lookup radius adjusted to %.1e",
                       pmap -> maxDist0);
               error(WARNING, errmsg);
#endif
            }
#ifdef PMAP_LOOKUP_REDO
            else {
               sprintf(errmsg, "max photon lookup radius clamped to %.1e",
                       pmap -> maxDist0);
               error(WARNING, errmsg);
            }
#endif
         }
         
         /* Reset successful lookup counter */
         pmap -> numLookups = 0;
      }
      else {
         /* Bail out after warning if maxDist is fixed */
         if (maxDistFix > 0)
            return;
            
         /* Increment successful lookup counter and reduce max search radius if
          * wraparound */
         pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT;
         if (!pmap -> numLookups)
            pmap -> maxDist0 *= PMAP_MAXDIST_DEC;
            
         redo = 0;
      }
      
   } while (redo);
}


void find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
{
   pmap -> maxDist2 = thescene.cusize;  /* ? */
   
#ifdef PMAP_OOC
   OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
#else
   kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
#endif                  
}


void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon)
{
#ifdef PMAP_OOC
   if (OOC_GetPhoton(pmap, idx, photon))      
#else
   if (kdT_GetPhoton(pmap, idx, photon))
#endif
      error(INTERNAL, "failed photon lookup");
}


Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
{
#ifdef PMAP_OOC
   return OOC_GetNearestPhoton(squeue, idx);
#else
   return kdT_GetNearestPhoton(squeue, idx);
#endif
}


PhotonIdx firstPhoton (const PhotonMap *pmap)
{
#ifdef PMAP_OOC
   return OOC_FirstPhoton(pmap);
#else
   return kdT_FirstPhoton(pmap);
#endif
}


void deletePhotons (PhotonMap* pmap)
{
#ifdef PMAP_OOC
   OOC_Delete(&pmap -> store);
#else
   kdT_Delete(&pmap -> store);
#endif

   free(pmap -> squeue.node);
   free(pmap -> biasCompHist);
   
   pmap -> numPhotons = pmap -> minGather = pmap -> maxGather =
      pmap -> squeue.len = pmap -> squeue.tail = 0;
}
Revision:	2.18
Committed:	Mon Aug 14 21:12:10 2017 UTC (6 years, 9 months ago) by rschregle
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R1
Changes since 2.17:	+97 -74 lines
Log Message:	Updated photon map code for Windows; no multproc or ooC for now
#	User	Rev	Content
1	greg	2.11	#ifndef lint
2	rschregle	2.18	static const char RCSid[] = "$Id: pmapdata.c,v 2.17 2016/09/29 20:51:07 greg Exp $";
3	greg	2.11	#endif
4	rschregle	2.15
5	greg	2.1	/*
6	rschregle	2.15	=========================================================================
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	greg	2.1
14			Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
15			(c) Fraunhofer Institute for Solar Energy Systems,
16	rschregle	2.3	(c) Lucerne University of Applied Sciences and Arts,
17	rschregle	2.15	supported by the Swiss National Science Foundation (SNSF, #147053)
18			==========================================================================
19	greg	2.1
20	rschregle	2.18	$Id: pmapdata.c,v 2.17 2016/09/29 20:51:07 greg Exp $
21	greg	2.1	*/
22
23
24
25	rschregle	2.18	#include "pmapdata.h"
26	greg	2.1	#include "pmaprand.h"
27			#include "pmapmat.h"
28			#include "otypes.h"
29			#include "source.h"
30			#include "rcontrib.h"
31	rschregle	2.7	#include "random.h"
32	greg	2.1
33
34
35			PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
36			NULL, NULL, NULL, NULL, NULL, NULL
37			};
38
39
40
41	rschregle	2.15	/* 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	greg	2.1	void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
51			/* Init photon map 'n' stuff... */
52			{
53			if (!pmap)
54			return;
55
56	rschregle	2.15	pmap -> numPhotons = 0;
57	greg	2.1	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;
60			pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0;
61			pmap -> gatherTolerance = gatherTolerance;
62			pmap -> minError = pmap -> maxError = pmap -> rmsError = 0;
63			pmap -> numDensity = 0;
64			pmap -> distribRatio = 1;
65			pmap -> type = t;
66	rschregle	2.15	pmap -> squeue.node = NULL;
67			pmap -> squeue.len = 0;
68	greg	2.1
69			/* Init local RNG state */
70			pmap -> randState [0] = 10243;
71			pmap -> randState [1] = 39829;
72			pmap -> randState [2] = 9433;
73			pmapSeed(randSeed, pmap -> randState);
74
75			/* Set up type-specific photon lookup callback */
76			pmap -> lookup = pmapLookup [t];
77
78	rschregle	2.15	/* Mark primary photon ray as unused */
79			pmap -> lastPrimary.srcIdx = -1;
80			pmap -> numPrimary = 0;
81			pmap -> primaries = NULL;
82
83			/* Init storage */
84			pmap -> heap = NULL;
85			pmap -> heapBuf = NULL;
86			pmap -> heapBufLen = 0;
87			#ifdef PMAP_OOC
88			OOC_Null(&pmap -> store);
89			#else
90			kdT_Null(&pmap -> store);
91			#endif
92	greg	2.1	}
93
94
95
96	rschregle	2.15	void initPhotonHeap (PhotonMap *pmap)
97	greg	2.1	{
98	rschregle	2.15	int fdFlags;
99	greg	2.1
100	rschregle	2.15	if (!pmap)
101			error(INTERNAL, "undefined photon map in initPhotonHeap");
102	greg	2.1
103	rschregle	2.15	if (!pmap -> heap) {
104			/* Open heap file */
105	rschregle	2.18	mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME));
106			if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b")))
107	rschregle	2.15	error(SYSTEM, "failed opening heap file in initPhotonHeap");
108	rschregle	2.18
109	greg	2.16	#ifdef F_SETFL /* XXX is there an alternate needed for Windows? */
110	rschregle	2.15	fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
111			fcntl(fileno(pmap -> heap), F_SETFL, fdFlags \| O_APPEND);
112	rschregle	2.18	#endif/* ftruncate(fileno(pmap -> heap), 0); */
113	greg	2.1	}
114	rschregle	2.15	}
115
116
117
118			void flushPhotonHeap (PhotonMap *pmap)
119			{
120			int fd;
121			const unsigned long len = pmap -> heapBufLen * sizeof(Photon);
122
123			if (!pmap)
124			error(INTERNAL, "undefined photon map in flushPhotonHeap");
125
126	rschregle	2.18	if (!pmap -> heap \|\| !pmap -> heapBuf) {
127			/* Silently ignore undefined heap
128			error(INTERNAL, "undefined heap in flushPhotonHeap"); */
129			return;
130			}
131	rschregle	2.15
132			/* Atomically seek and write block to heap */
133			/* !!! Unbuffered I/O via pwrite() avoids potential race conditions
134			* !!! and buffer corruption which can occur with lseek()/fseek()
135			* !!! followed by write()/fwrite(). */
136			fd = fileno(pmap -> heap);
137	greg	2.1
138	rschregle	2.15	#ifdef DEBUG_PMAP
139			sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(),
140			pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon));
141			eputs(errmsg);
142			#endif
143
144			/if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) /
145			if (write(fd, pmap -> heapBuf, len) != len)
146			error(SYSTEM, "failed append to heap file in flushPhotonHeap");
147	greg	2.17
148	rschregle	2.18	#if NIX
149	rschregle	2.15	if (fsync(fd))
150			error(SYSTEM, "failed fsync in flushPhotonHeap");
151	greg	2.17	#endif
152	rschregle	2.18
153	rschregle	2.15	pmap -> heapBufLen = 0;
154			}
155	greg	2.1
156	rschregle	2.15
157
158	rschregle	2.18	#ifdef DEBUG_PMAP
159	rschregle	2.15	static int checkPhotonHeap (FILE *file)
160			/* Check heap for nonsensical or duplicate photons */
161			{
162			Photon p, lastp;
163			int i, dup;
164
165			rewind(file);
166			memset(&lastp, 0, sizeof(lastp));
167	greg	2.1
168	rschregle	2.15	while (fread(&p, sizeof(p), 1, file)) {
169			dup = 1;
170
171			for (i = 0; i <= 2; i++) {
172			if (p.pos [i] < thescene.cuorg [i] \|\|
173			p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
174
175			sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
176			p.pos [0], p.pos [1], p.pos [2]);
177			error(WARNING, errmsg);
178			}
179
180			dup &= p.pos [i] == lastp.pos [i];
181			}
182
183			if (dup) {
184			sprintf(errmsg,
185			"consecutive duplicate photon in heap at [%f, %f, %f]\n",
186			p.pos [0], p.pos [1], p.pos [2]);
187			error(WARNING, errmsg);
188			}
189			}
190
191			return 0;
192	greg	2.1	}
193	rschregle	2.15	#endif
194	greg	2.1
195
196
197	rschregle	2.15	int newPhoton (PhotonMap* pmap, const RAY* ray)
198	greg	2.1	{
199	rschregle	2.18	unsigned i, inROI = 0;
200	rschregle	2.15	Photon photon;
201	greg	2.1	COLOR photonFlux;
202
203			/* Account for distribution ratio */
204			if (!pmap \|\| pmapRandom(pmap -> randState) > pmap -> distribRatio)
205	rschregle	2.15	return -1;
206	greg	2.1
207			/* Don't store on sources */
208			if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
209	rschregle	2.15	return -1;
210	greg	2.1
211	rschregle	2.18	/* Store photon if within a region of interest (for ze Ecksperts!) */
212			if (!pmapNumROI \|\| !pmapROI)
213			inROI = 1;
214			else {
215			for (i = 0; !inROI && i < pmapNumROI; i++)
216			inROI = (ray -> rop [0] >= pmapROI [i].min [0] &&
217			ray -> rop [0] <= pmapROI [i].max [0] &&
218			ray -> rop [1] >= pmapROI [i].min [1] &&
219			ray -> rop [1] <= pmapROI [i].max [1] &&
220			ray -> rop [2] >= pmapROI [i].min [2] &&
221			ray -> rop [2] <= pmapROI [i].max [2]);
222			}
223
224			if (inROI) {
225	greg	2.1	/* Adjust flux according to distribution ratio and ray weight */
226			copycolor(photonFlux, ray -> rcol);
227			scalecolor(photonFlux,
228			ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
229			: 1));
230	rschregle	2.15	setPhotonFlux(&photon, photonFlux);
231	greg	2.1
232			/* Set photon position and flags */
233	rschregle	2.15	VCOPY(photon.pos, ray -> rop);
234			photon.flags = 0;
235			photon.caustic = PMAP_CAUSTICRAY(ray);
236
237			/* Set contrib photon's primary ray and subprocess index (the latter
238			* to linearise the primary ray indices after photon distribution is
239			* complete). Also set primary ray's source index, thereby marking it
240			* as used. */
241			if (isContribPmap(pmap)) {
242			photon.primary = pmap -> numPrimary;
243			photon.proc = PMAP_GETRAYPROC(ray);
244			pmap -> lastPrimary.srcIdx = ray -> rsrc;
245			}
246			else photon.primary = 0;
247
248			/* Set normal */
249			for (i = 0; i <= 2; i++)
250			photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
251			: ray -> ron [i]);
252
253			if (!pmap -> heapBuf) {
254			/* Lazily allocate heap buffa */
255	rschregle	2.18	#if NIX
256			/* Randomise buffa size to temporally decorellate flushes in
257			* multiprocessing mode */
258	rschregle	2.15	srandom(randSeed + getpid());
259			pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
260			#else
261	rschregle	2.18	/* Randomisation disabled for single processes on WIN; also useful
262			* for reproducability during debugging */
263	rschregle	2.15	pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
264			#endif
265			if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
266			error(SYSTEM, "failed heap buffer allocation in newPhoton");
267			pmap -> heapBufLen = 0;
268			}
269
270			/* Photon initialised; now append to heap buffa */
271			memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
272
273			if (++pmap -> heapBufLen >= pmap -> heapBufSize)
274			/* Heap buffa full, flush to heap file */
275			flushPhotonHeap(pmap);
276	greg	2.1
277	rschregle	2.15	pmap -> numPhotons++;
278	greg	2.1	}
279
280	rschregle	2.15	return 0;
281	greg	2.1	}
282
283
284
285	rschregle	2.15	void buildPhotonMap (PhotonMap pmap, double photonFlux,
286			PhotonPrimaryIdx *primaryOfs, unsigned nproc)
287			{
288			unsigned long n, nCheck = 0;
289			unsigned i;
290			Photon *p;
291			COLOR flux;
292	rschregle	2.18	char nuHeapFname [sizeof(PMAP_TMPFNAME)];
293	rschregle	2.15	FILE *nuHeap;
294			/* Need double here to reduce summation errors */
295			double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
296			FVECT d;
297
298			if (!pmap)
299			error(INTERNAL, "undefined photon map in buildPhotonMap");
300
301			/* Get number of photons from heapfile size */
302	rschregle	2.18	if (fseek(pmap -> heap, 0, SEEK_END) < 0)
303			error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap");
304	rschregle	2.15	pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);
305
306			if (!pmap -> numPhotons)
307			error(INTERNAL, "empty photon map in buildPhotonMap");
308
309			if (!pmap -> heap)
310			error(INTERNAL, "no heap in buildPhotonMap");
311
312			#ifdef DEBUG_PMAP
313			eputs("Checking photon heap consistency...\n");
314			checkPhotonHeap(pmap -> heap);
315
316			sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons);
317			eputs(errmsg);
318			#endif
319	rschregle	2.18
320	rschregle	2.15	/* Allocate heap buffa */
321			if (!pmap -> heapBuf) {
322			pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
323			pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon));
324			if (!pmap -> heapBuf)
325			error(SYSTEM, "failed to allocate postprocessed photon heap in"
326			"buildPhotonMap");
327	greg	2.1	}
328
329	rschregle	2.15	/* We REALLY don't need yet another @%&*! heap just to hold the scaled
330			* photons, but can't think of a quicker fix... */
331	rschregle	2.18	mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME));
332			if (!(nuHeap = fopen(nuHeapFname, "w+b")))
333	rschregle	2.15	error(SYSTEM, "failed to open postprocessed photon heap in "
334			"buildPhotonMap");
335
336			rewind(pmap -> heap);
337
338			#ifdef DEBUG_PMAP
339			eputs("Postprocessing photons...\n");
340			#endif
341	greg	2.1
342	rschregle	2.18	while (!feof(pmap -> heap)) {
343			#ifdef DEBUG_PMAP
344			printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap));
345			#endif
346	rschregle	2.15	pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
347	rschregle	2.18	pmap -> heapBufSize, pmap -> heap);
348			#ifdef DEBUG_PMAP
349			printf("Got %lu\n", pmap -> heapBufLen);
350			#endif
351
352			if (ferror(pmap -> heap))
353			error(SYSTEM, "failed to read photon heap in buildPhotonMap");
354	rschregle	2.15
355	rschregle	2.18	for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
356			/* Update min and max pos and set photon flux */
357			for (i = 0; i <= 2; i++) {
358			if (p -> pos [i] < pmap -> minPos [i])
359			pmap -> minPos [i] = p -> pos [i];
360			else if (p -> pos [i] > pmap -> maxPos [i])
361			pmap -> maxPos [i] = p -> pos [i];
362
363			/* Update centre of gravity with photon position */
364			CoG [i] += p -> pos [i];
365			}
366
367			if (primaryOfs)
368			/* Linearise photon primary index from subprocess index using the
369			* per-subprocess offsets in primaryOfs */
370			p -> primary += primaryOfs [p -> proc];
371
372			/* Scale photon's flux (hitherto normalised to 1 over RGB); in
373			* case of a contrib photon map, this is done per light source,
374			* and photonFlux is assumed to be an array */
375			getPhotonFlux(p, flux);
376
377			if (photonFlux) {
378			scalecolor(flux, photonFlux [isContribPmap(pmap) ?
379			photonSrcIdx(pmap, p) : 0]);
380			setPhotonFlux(p, flux);
381	greg	2.1	}
382	rschregle	2.18
383			/* Update average photon flux; need a double here */
384			addcolor(avgFlux, flux);
385			}
386	greg	2.1
387	rschregle	2.18	/* Write modified photons to new heap */
388			fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
389	rschregle	2.15
390	rschregle	2.18	if (ferror(nuHeap))
391			error(SYSTEM, "failed postprocessing photon flux in "
392			"buildPhotonMap");
393	rschregle	2.15
394			nCheck += pmap -> heapBufLen;
395			}
396
397			#ifdef DEBUG_PMAP
398			if (nCheck < pmap -> numPhotons)
399			error(INTERNAL, "truncated photon heap in buildPhotonMap");
400			#endif
401
402			/* Finalise average photon flux */
403			scalecolor(avgFlux, 1.0 / pmap -> numPhotons);
404			copycolor(pmap -> photonFlux, avgFlux);
405
406			/* Average photon positions to get centre of gravity */
407			for (i = 0; i < 3; i++)
408			pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons;
409
410			rewind(pmap -> heap);
411
412			/* Compute average photon distance to centre of gravity */
413			while (!feof(pmap -> heap)) {
414			pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
415	rschregle	2.18	pmap -> heapBufSize, pmap -> heap);
416	rschregle	2.15
417	rschregle	2.18	for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
418			VSUB(d, p -> pos, CoG);
419			CoGdist += DOT(d, d);
420			}
421	rschregle	2.15	}
422
423			pmap -> CoGdist = CoGdist /= pmap -> numPhotons;
424
425	rschregle	2.18	/* Swap heaps, discarding unscaled photons */
426	rschregle	2.15	fclose(pmap -> heap);
427	rschregle	2.18	unlink(pmap -> heapFname);
428	rschregle	2.15	pmap -> heap = nuHeap;
429	rschregle	2.18	strcpy(pmap -> heapFname, nuHeapFname);
430	rschregle	2.15
431			#ifdef PMAP_OOC
432			OOC_BuildPhotonMap(pmap, nproc);
433			#else
434			kdT_BuildPhotonMap(pmap);
435			#endif
436
437			/* Trash heap and its buffa */
438			free(pmap -> heapBuf);
439			fclose(pmap -> heap);
440	rschregle	2.18	unlink(pmap -> heapFname);
441	rschregle	2.15	pmap -> heap = NULL;
442			pmap -> heapBuf = NULL;
443	greg	2.1	}
444
445
446
447			/* Dynamic max photon search radius increase and reduction factors */
448			#define PMAP_MAXDIST_INC 4
449			#define PMAP_MAXDIST_DEC 0.9
450
451			/* Num successful lookups before reducing in max search radius */
452			#define PMAP_MAXDIST_CNT 1000
453
454			/* Threshold below which we assume increasing max radius won't help */
455			#define PMAP_SHORT_LOOKUP_THRESH 1
456
457	rschregle	2.8	/* Coefficient for adaptive maximum search radius */
458			#define PMAP_MAXDIST_COEFF 100
459
460	greg	2.1	void findPhotons (PhotonMap* pmap, const RAY* ray)
461			{
462			int redo = 0;
463
464	rschregle	2.15	if (!pmap -> squeue.len) {
465	greg	2.1	/* Lazy init priority queue */
466	rschregle	2.15	#ifdef PMAP_OOC
467			OOC_InitFindPhotons(pmap);
468			#else
469			kdT_InitFindPhotons(pmap);
470			#endif
471	greg	2.1	pmap -> minGathered = pmap -> maxGather;
472			pmap -> maxGathered = pmap -> minGather;
473			pmap -> totalGathered = 0;
474			pmap -> numLookups = pmap -> numShortLookups = 0;
475			pmap -> shortLookupPct = 0;
476			pmap -> minError = FHUGE;
477			pmap -> maxError = -FHUGE;
478			pmap -> rmsError = 0;
479	rschregle	2.9	/* SQUARED max search radius limit is based on avg photon distance to
480	rschregle	2.8	* centre of gravity, unless fixed by user (maxDistFix > 0) */
481	rschregle	2.15	pmap -> maxDist0 = pmap -> maxDist2Limit =
482	rschregle	2.9	maxDistFix > 0 ? maxDistFix * maxDistFix
483	rschregle	2.15	: PMAP_MAXDIST_COEFF * pmap -> squeue.len *
484			pmap -> CoGdist / pmap -> numPhotons;
485	greg	2.1	}
486
487			do {
488	rschregle	2.15	pmap -> squeue.tail = 0;
489			pmap -> maxDist2 = pmap -> maxDist0;
490	greg	2.1
491			/* Search position is ray -> rorg for volume photons, since we have no
492			intersection point. Normals are ignored -- these are incident
493			directions). */
494			if (isVolumePmap(pmap)) {
495	rschregle	2.15	#ifdef PMAP_OOC
496			OOC_FindPhotons(pmap, ray -> rorg, NULL);
497			#else
498			kdT_FindPhotons(pmap, ray -> rorg, NULL);
499			#endif
500	greg	2.1	}
501			else {
502	rschregle	2.15	#ifdef PMAP_OOC
503			OOC_FindPhotons(pmap, ray -> rop, ray -> ron);
504			#else
505			kdT_FindPhotons(pmap, ray -> rop, ray -> ron);
506			#endif
507	greg	2.1	}
508	rschregle	2.3
509	rschregle	2.15	#ifdef PMAP_LOOKUP_INFO
510			fprintf(stderr, "%d/%d %s photons found within radius %.3f "
511			"at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail,
512			pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2),
513			ray -> rop [0], ray -> rop [1], ray -> rop [2],
514			ray -> ro ? ray -> ro -> oname : "<null>");
515			#endif
516
517			if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) {
518	greg	2.1	/* Short lookup; too few photons found */
519	rschregle	2.15	if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) {
520	greg	2.1	/* Ignore short lookups which return fewer than
521			* PMAP_SHORT_LOOKUP_THRESH photons under the assumption there
522			* really are no photons in the vicinity, and increasing the max
523			* search radius therefore won't help */
524	rschregle	2.8	#ifdef PMAP_LOOKUP_WARN
525	greg	2.1	sprintf(errmsg,
526			"%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
527	rschregle	2.15	pmap -> squeue.tail, pmap -> squeue.len,
528			pmapName [pmap -> type],
529			ray -> rop [0], ray -> rop [1], ray -> rop [2],
530	greg	2.1	ray -> ro ? ray -> ro -> oname : "<null>");
531			error(WARNING, errmsg);
532	rschregle	2.8	#endif
533	rschregle	2.3
534	rschregle	2.8	/* Bail out after warning if maxDist is fixed */
535			if (maxDistFix > 0)
536			return;
537
538	rschregle	2.15	if (pmap -> maxDist0 < pmap -> maxDist2Limit) {
539	greg	2.1	/* Increase max search radius if below limit & redo search */
540			pmap -> maxDist0 *= PMAP_MAXDIST_INC;
541	rschregle	2.8	#ifdef PMAP_LOOKUP_REDO
542	greg	2.1	redo = 1;
543	rschregle	2.8	#endif
544			#ifdef PMAP_LOOKUP_WARN
545	greg	2.1	sprintf(errmsg,
546			redo ? "restarting photon lookup with max radius %.1e"
547			: "max photon lookup radius adjusted to %.1e",
548	rschregle	2.15	pmap -> maxDist0);
549	greg	2.1	error(WARNING, errmsg);
550	rschregle	2.8	#endif
551	greg	2.1	}
552	rschregle	2.8	#ifdef PMAP_LOOKUP_REDO
553	greg	2.1	else {
554			sprintf(errmsg, "max photon lookup radius clamped to %.1e",
555	rschregle	2.15	pmap -> maxDist0);
556	greg	2.1	error(WARNING, errmsg);
557			}
558	rschregle	2.8	#endif
559	greg	2.1	}
560
561			/* Reset successful lookup counter */
562			pmap -> numLookups = 0;
563	rschregle	2.3	}
564	greg	2.1	else {
565	rschregle	2.8	/* Bail out after warning if maxDist is fixed */
566			if (maxDistFix > 0)
567			return;
568
569	greg	2.1	/* Increment successful lookup counter and reduce max search radius if
570			* wraparound */
571			pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT;
572			if (!pmap -> numLookups)
573			pmap -> maxDist0 *= PMAP_MAXDIST_DEC;
574
575			redo = 0;
576			}
577	rschregle	2.8
578	greg	2.1	} while (redo);
579			}
580
581
582
583	rschregle	2.15	void find1Photon (PhotonMap pmap, const RAY ray, Photon *photon)
584			{
585			pmap -> maxDist2 = thescene.cusize; /* ? */
586	greg	2.1
587	rschregle	2.15	#ifdef PMAP_OOC
588			OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
589			#else
590			kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
591			#endif
592	greg	2.1	}
593
594
595
596	rschregle	2.15	void getPhoton (PhotonMap pmap, PhotonIdx idx, Photon photon)
597	greg	2.1	{
598	rschregle	2.15	#ifdef PMAP_OOC
599	rschregle	2.18	if (OOC_GetPhoton(pmap, idx, photon))
600	rschregle	2.15	#else
601			if (kdT_GetPhoton(pmap, idx, photon))
602			#endif
603			error(INTERNAL, "failed photon lookup");
604	greg	2.1	}
605
606
607
608	rschregle	2.15	Photon getNearestPhoton (const PhotonSearchQueue squeue, PhotonIdx idx)
609			{
610			#ifdef PMAP_OOC
611			return OOC_GetNearestPhoton(squeue, idx);
612			#else
613			return kdT_GetNearestPhoton(squeue, idx);
614			#endif
615	greg	2.1	}
616
617
618
619	rschregle	2.15	PhotonIdx firstPhoton (const PhotonMap *pmap)
620			{
621			#ifdef PMAP_OOC
622			return OOC_FirstPhoton(pmap);
623			#else
624			return kdT_FirstPhoton(pmap);
625			#endif
626	greg	2.1	}
627
628
629
630			void deletePhotons (PhotonMap* pmap)
631			{
632	rschregle	2.15	#ifdef PMAP_OOC
633			OOC_Delete(&pmap -> store);
634			#else
635			kdT_Delete(&pmap -> store);
636			#endif
637
638			free(pmap -> squeue.node);
639	greg	2.1	free(pmap -> biasCompHist);
640
641	rschregle	2.15	pmap -> numPhotons = pmap -> minGather = pmap -> maxGather =
642			pmap -> squeue.len = pmap -> squeue.tail = 0;
643	greg	2.1	}