src/rt/pmapdata.c

#ifndef lint
static const char RCSid[] = "$Id: pmapdata.c,v 2.18 2017/08/14 21:12:10 rschregle 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.18 2017/08/14 21:12:10 rschregle 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

/* Ambient include/exclude set (from ambient.c) */
#ifndef  MAXASET
   #define MAXASET  4095
#endif
extern OBJECT ambset [MAXASET+1];


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;
   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;

   /*  if modifier in include/exclude set */
   if (ambincl != -1 && ray -> ro && 
       ambincl != inset(ambset, ray -> ro -> omod))
      return -1;

   if (pmapNumROI && pmapROI) {      
      unsigned inROI = 0;
      
      /* Store photon if within a region of interest (for ze Ecksperts!) */
      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)
         return -1;
   }
    
   /* 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.19
Committed:	Tue Mar 20 19:55:33 2018 UTC (6 years, 1 month ago) by rschregle
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2
Changes since 2.18:	+68 -57 lines
Log Message:	Added -ae/-ai ambient exclude options to mkpmap, cleaned up opt parsing.
#	User	Rev	Content
1	greg	2.11	#ifndef lint
2	rschregle	2.19	static const char RCSid[] = "$Id: pmapdata.c,v 2.18 2017/08/14 21:12:10 rschregle 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.19	$Id: pmapdata.c,v 2.18 2017/08/14 21:12:10 rschregle 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	rschregle	2.19	/* Ambient include/exclude set (from ambient.c) */
49			#ifndef MAXASET
50			#define MAXASET 4095
51			#endif
52			extern OBJECT ambset [MAXASET+1];
53
54	rschregle	2.15
55
56	greg	2.1	void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
57			/* Init photon map 'n' stuff... */
58			{
59			if (!pmap)
60			return;
61
62	rschregle	2.15	pmap -> numPhotons = 0;
63	greg	2.1	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;
66			pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0;
67			pmap -> gatherTolerance = gatherTolerance;
68			pmap -> minError = pmap -> maxError = pmap -> rmsError = 0;
69			pmap -> numDensity = 0;
70			pmap -> distribRatio = 1;
71			pmap -> type = t;
72	rschregle	2.15	pmap -> squeue.node = NULL;
73			pmap -> squeue.len = 0;
74	greg	2.1
75			/* Init local RNG state */
76			pmap -> randState [0] = 10243;
77			pmap -> randState [1] = 39829;
78			pmap -> randState [2] = 9433;
79			pmapSeed(randSeed, pmap -> randState);
80
81			/* Set up type-specific photon lookup callback */
82			pmap -> lookup = pmapLookup [t];
83
84	rschregle	2.15	/* 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	greg	2.1	}
99
100
101
102	rschregle	2.15	void initPhotonHeap (PhotonMap *pmap)
103	greg	2.1	{
104	rschregle	2.15	int fdFlags;
105	greg	2.1
106	rschregle	2.15	if (!pmap)
107			error(INTERNAL, "undefined photon map in initPhotonHeap");
108	greg	2.1
109	rschregle	2.15	if (!pmap -> heap) {
110			/* Open heap file */
111	rschregle	2.18	mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME));
112			if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b")))
113	rschregle	2.15	error(SYSTEM, "failed opening heap file in initPhotonHeap");
114	rschregle	2.18
115	greg	2.16	#ifdef F_SETFL /* XXX is there an alternate needed for Windows? */
116	rschregle	2.15	fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
117			fcntl(fileno(pmap -> heap), F_SETFL, fdFlags \| O_APPEND);
118	rschregle	2.18	#endif/* ftruncate(fileno(pmap -> heap), 0); */
119	greg	2.1	}
120	rschregle	2.15	}
121
122
123
124			void flushPhotonHeap (PhotonMap *pmap)
125			{
126			int fd;
127			const unsigned long len = pmap -> heapBufLen * sizeof(Photon);
128
129			if (!pmap)
130			error(INTERNAL, "undefined photon map in flushPhotonHeap");
131
132	rschregle	2.18	if (!pmap -> heap \|\| !pmap -> heapBuf) {
133			/* Silently ignore undefined heap
134			error(INTERNAL, "undefined heap in flushPhotonHeap"); */
135			return;
136			}
137	rschregle	2.15
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	greg	2.1
144	rschregle	2.15	#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			/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	greg	2.17
154	rschregle	2.18	#if NIX
155	rschregle	2.15	if (fsync(fd))
156			error(SYSTEM, "failed fsync in flushPhotonHeap");
157	greg	2.17	#endif
158	rschregle	2.18
159	rschregle	2.15	pmap -> heapBufLen = 0;
160			}
161	greg	2.1
162	rschregle	2.15
163
164	rschregle	2.18	#ifdef DEBUG_PMAP
165	rschregle	2.15	static int checkPhotonHeap (FILE *file)
166			/* Check heap for nonsensical or duplicate photons */
167			{
168			Photon p, lastp;
169			int i, dup;
170
171			rewind(file);
172			memset(&lastp, 0, sizeof(lastp));
173	greg	2.1
174	rschregle	2.15	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	greg	2.1	}
199	rschregle	2.15	#endif
200	greg	2.1
201
202
203	rschregle	2.15	int newPhoton (PhotonMap* pmap, const RAY* ray)
204	greg	2.1	{
205	rschregle	2.19	unsigned i;
206	rschregle	2.15	Photon photon;
207	greg	2.1	COLOR photonFlux;
208
209			/* Account for distribution ratio */
210			if (!pmap \|\| pmapRandom(pmap -> randState) > pmap -> distribRatio)
211	rschregle	2.15	return -1;
212	greg	2.1
213			/* Don't store on sources */
214			if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
215	rschregle	2.15	return -1;
216	greg	2.1
217	rschregle	2.19	/* if modifier in include/exclude 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			/* Store photon if within a region of interest (for ze Ecksperts!) */
226	rschregle	2.18	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	rschregle	2.19	if (!inROI)
234			return -1;
235	rschregle	2.18	}
236
237	rschregle	2.19	/* 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			/* Set photon position and flags */
245			VCOPY(photon.pos, ray -> rop);
246			photon.flags = 0;
247			photon.caustic = PMAP_CAUSTICRAY(ray);
248
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	rschregle	2.15
265	rschregle	2.19	if (!pmap -> heapBuf) {
266			/* Lazily allocate heap buffa */
267	rschregle	2.18	#if NIX
268	rschregle	2.19	/* 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	rschregle	2.15	#endif
277	rschregle	2.19	if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
278			error(SYSTEM, "failed heap buffer allocation in newPhoton");
279			pmap -> heapBufLen = 0;
280			}
281	rschregle	2.15
282	rschregle	2.19	/* 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	greg	2.1
289	rschregle	2.19	pmap -> numPhotons++;
290	greg	2.1
291	rschregle	2.15	return 0;
292	greg	2.1	}
293
294
295
296	rschregle	2.15	void buildPhotonMap (PhotonMap pmap, double photonFlux,
297			PhotonPrimaryIdx *primaryOfs, unsigned nproc)
298			{
299			unsigned long n, nCheck = 0;
300			unsigned i;
301			Photon *p;
302			COLOR flux;
303	rschregle	2.18	char nuHeapFname [sizeof(PMAP_TMPFNAME)];
304	rschregle	2.15	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			if (!pmap)
310			error(INTERNAL, "undefined photon map in buildPhotonMap");
311
312			/* Get number of photons from heapfile size */
313	rschregle	2.18	if (fseek(pmap -> heap, 0, SEEK_END) < 0)
314			error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap");
315	rschregle	2.15	pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);
316
317			if (!pmap -> numPhotons)
318			error(INTERNAL, "empty photon map in buildPhotonMap");
319
320			if (!pmap -> heap)
321			error(INTERNAL, "no heap in buildPhotonMap");
322
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	rschregle	2.18
331	rschregle	2.15	/* 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	greg	2.1	}
339
340	rschregle	2.15	/* We REALLY don't need yet another @%&*! heap just to hold the scaled
341			* photons, but can't think of a quicker fix... */
342	rschregle	2.18	mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME));
343			if (!(nuHeap = fopen(nuHeapFname, "w+b")))
344	rschregle	2.15	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	greg	2.1
353	rschregle	2.18	while (!feof(pmap -> heap)) {
354			#ifdef DEBUG_PMAP
355			printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap));
356			#endif
357	rschregle	2.15	pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
358	rschregle	2.18	pmap -> heapBufSize, pmap -> heap);
359			#ifdef DEBUG_PMAP
360			printf("Got %lu\n", pmap -> heapBufLen);
361			#endif
362
363			if (ferror(pmap -> heap))
364			error(SYSTEM, "failed to read photon heap in buildPhotonMap");
365	rschregle	2.15
366	rschregle	2.18	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			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			/* 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	greg	2.1	}
393	rschregle	2.18
394			/* Update average photon flux; need a double here */
395			addcolor(avgFlux, flux);
396			}
397	greg	2.1
398	rschregle	2.18	/* Write modified photons to new heap */
399			fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
400	rschregle	2.15
401	rschregle	2.18	if (ferror(nuHeap))
402			error(SYSTEM, "failed postprocessing photon flux in "
403			"buildPhotonMap");
404	rschregle	2.15
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	rschregle	2.18	pmap -> heapBufSize, pmap -> heap);
427	rschregle	2.15
428	rschregle	2.18	for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
429			VSUB(d, p -> pos, CoG);
430			CoGdist += DOT(d, d);
431			}
432	rschregle	2.15	}
433
434			pmap -> CoGdist = CoGdist /= pmap -> numPhotons;
435
436	rschregle	2.18	/* Swap heaps, discarding unscaled photons */
437	rschregle	2.15	fclose(pmap -> heap);
438	rschregle	2.18	unlink(pmap -> heapFname);
439	rschregle	2.15	pmap -> heap = nuHeap;
440	rschregle	2.18	strcpy(pmap -> heapFname, nuHeapFname);
441	rschregle	2.15
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	rschregle	2.18	unlink(pmap -> heapFname);
452	rschregle	2.15	pmap -> heap = NULL;
453			pmap -> heapBuf = NULL;
454	greg	2.1	}
455
456
457
458			/* Dynamic max photon search radius increase and reduction factors */
459			#define PMAP_MAXDIST_INC 4
460			#define PMAP_MAXDIST_DEC 0.9
461
462			/* Num successful lookups before reducing in max search radius */
463			#define PMAP_MAXDIST_CNT 1000
464
465			/* Threshold below which we assume increasing max radius won't help */
466			#define PMAP_SHORT_LOOKUP_THRESH 1
467
468	rschregle	2.8	/* Coefficient for adaptive maximum search radius */
469			#define PMAP_MAXDIST_COEFF 100
470
471	greg	2.1	void findPhotons (PhotonMap* pmap, const RAY* ray)
472			{
473			int redo = 0;
474
475	rschregle	2.15	if (!pmap -> squeue.len) {
476	greg	2.1	/* Lazy init priority queue */
477	rschregle	2.15	#ifdef PMAP_OOC
478			OOC_InitFindPhotons(pmap);
479			#else
480			kdT_InitFindPhotons(pmap);
481			#endif
482	greg	2.1	pmap -> minGathered = pmap -> maxGather;
483			pmap -> maxGathered = pmap -> minGather;
484			pmap -> totalGathered = 0;
485			pmap -> numLookups = pmap -> numShortLookups = 0;
486			pmap -> shortLookupPct = 0;
487			pmap -> minError = FHUGE;
488			pmap -> maxError = -FHUGE;
489			pmap -> rmsError = 0;
490	rschregle	2.9	/* SQUARED max search radius limit is based on avg photon distance to
491	rschregle	2.8	* centre of gravity, unless fixed by user (maxDistFix > 0) */
492	rschregle	2.15	pmap -> maxDist0 = pmap -> maxDist2Limit =
493	rschregle	2.9	maxDistFix > 0 ? maxDistFix * maxDistFix
494	rschregle	2.15	: PMAP_MAXDIST_COEFF * pmap -> squeue.len *
495			pmap -> CoGdist / pmap -> numPhotons;
496	greg	2.1	}
497
498			do {
499	rschregle	2.15	pmap -> squeue.tail = 0;
500			pmap -> maxDist2 = pmap -> maxDist0;
501	greg	2.1
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	rschregle	2.15	#ifdef PMAP_OOC
507			OOC_FindPhotons(pmap, ray -> rorg, NULL);
508			#else
509			kdT_FindPhotons(pmap, ray -> rorg, NULL);
510			#endif
511	greg	2.1	}
512			else {
513	rschregle	2.15	#ifdef PMAP_OOC
514			OOC_FindPhotons(pmap, ray -> rop, ray -> ron);
515			#else
516			kdT_FindPhotons(pmap, ray -> rop, ray -> ron);
517			#endif
518	greg	2.1	}
519	rschregle	2.3
520	rschregle	2.15	#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	greg	2.1	/* Short lookup; too few photons found */
530	rschregle	2.15	if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) {
531	greg	2.1	/* 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
534			* search radius therefore won't help */
535	rschregle	2.8	#ifdef PMAP_LOOKUP_WARN
536	greg	2.1	sprintf(errmsg,
537			"%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
538	rschregle	2.15	pmap -> squeue.tail, pmap -> squeue.len,
539			pmapName [pmap -> type],
540			ray -> rop [0], ray -> rop [1], ray -> rop [2],
541	greg	2.1	ray -> ro ? ray -> ro -> oname : "<null>");
542			error(WARNING, errmsg);
543	rschregle	2.8	#endif
544	rschregle	2.3
545	rschregle	2.8	/* Bail out after warning if maxDist is fixed */
546			if (maxDistFix > 0)
547			return;
548
549	rschregle	2.15	if (pmap -> maxDist0 < pmap -> maxDist2Limit) {
550	greg	2.1	/* Increase max search radius if below limit & redo search */
551			pmap -> maxDist0 *= PMAP_MAXDIST_INC;
552	rschregle	2.8	#ifdef PMAP_LOOKUP_REDO
553	greg	2.1	redo = 1;
554	rschregle	2.8	#endif
555			#ifdef PMAP_LOOKUP_WARN
556	greg	2.1	sprintf(errmsg,
557			redo ? "restarting photon lookup with max radius %.1e"
558			: "max photon lookup radius adjusted to %.1e",
559	rschregle	2.15	pmap -> maxDist0);
560	greg	2.1	error(WARNING, errmsg);
561	rschregle	2.8	#endif
562	greg	2.1	}
563	rschregle	2.8	#ifdef PMAP_LOOKUP_REDO
564	greg	2.1	else {
565			sprintf(errmsg, "max photon lookup radius clamped to %.1e",
566	rschregle	2.15	pmap -> maxDist0);
567	greg	2.1	error(WARNING, errmsg);
568			}
569	rschregle	2.8	#endif
570	greg	2.1	}
571
572			/* Reset successful lookup counter */
573			pmap -> numLookups = 0;
574	rschregle	2.3	}
575	greg	2.1	else {
576	rschregle	2.8	/* Bail out after warning if maxDist is fixed */
577			if (maxDistFix > 0)
578			return;
579
580	greg	2.1	/* Increment successful lookup counter and reduce max search radius if
581			* wraparound */
582			pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT;
583			if (!pmap -> numLookups)
584			pmap -> maxDist0 *= PMAP_MAXDIST_DEC;
585
586			redo = 0;
587			}
588	rschregle	2.8
589	greg	2.1	} while (redo);
590			}
591
592
593
594	rschregle	2.15	void find1Photon (PhotonMap pmap, const RAY ray, Photon *photon)
595			{
596			pmap -> maxDist2 = thescene.cusize; /* ? */
597	greg	2.1
598	rschregle	2.15	#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	greg	2.1	}
604
605
606
607	rschregle	2.15	void getPhoton (PhotonMap pmap, PhotonIdx idx, Photon photon)
608	greg	2.1	{
609	rschregle	2.15	#ifdef PMAP_OOC
610	rschregle	2.18	if (OOC_GetPhoton(pmap, idx, photon))
611	rschregle	2.15	#else
612			if (kdT_GetPhoton(pmap, idx, photon))
613			#endif
614			error(INTERNAL, "failed photon lookup");
615	greg	2.1	}
616
617
618
619	rschregle	2.15	Photon getNearestPhoton (const PhotonSearchQueue squeue, PhotonIdx idx)
620			{
621			#ifdef PMAP_OOC
622			return OOC_GetNearestPhoton(squeue, idx);
623			#else
624			return kdT_GetNearestPhoton(squeue, idx);
625			#endif
626	greg	2.1	}
627
628
629
630	rschregle	2.15	PhotonIdx firstPhoton (const PhotonMap *pmap)
631			{
632			#ifdef PMAP_OOC
633			return OOC_FirstPhoton(pmap);
634			#else
635			return kdT_FirstPhoton(pmap);
636			#endif
637	greg	2.1	}
638
639
640
641			void deletePhotons (PhotonMap* pmap)
642			{
643	rschregle	2.15	#ifdef PMAP_OOC
644			OOC_Delete(&pmap -> store);
645			#else
646			kdT_Delete(&pmap -> store);
647			#endif
648
649			free(pmap -> squeue.node);
650	greg	2.1	free(pmap -> biasCompHist);
651
652	rschregle	2.15	pmap -> numPhotons = pmap -> minGather = pmap -> maxGather =
653			pmap -> squeue.len = pmap -> squeue.tail = 0;
654	greg	2.1	}