src/rt/pmapooc.c

/* 
   ======================================================================
   Photon map interface to out-of-core octree

   Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
   (c) Lucerne University of Applied Sciences and Arts,
       supported by the Swiss National Science Foundation (SNSF, #147053)
   ======================================================================
   
   $Id: pmapooc.c,v 1.5 2018/11/08 00:54:07 greg Exp $
*/


#include "pmapdata.h"   /* Includes pmapooc.h */
#include "source.h"
#include "otspecial.h"
#include "oocsort.h"
#include "oocbuild.h"


/* Returns photon position as sorting key for OOC_Octree & friends (notably
 * for Morton code generation).
 * !!! Uses type conversion from float via TEMPORARY storage; 
 * !!! THIS IS NOT THREAD SAFE!
 * !!! RETURNED KEY PERSISTS ONLY IF COPIED BEFORE NEXT CALL! */
RREAL *OOC_PhotonKey (const void *p)
{
   static FVECT photonPos; /* Temp storage for type conversion */
   
   VCOPY(photonPos, ((Photon*)p) -> pos);
   return photonPos;
}
           

#ifdef DEBUG_OOC
static int OOC_CheckKeys (FILE *file, const OOC_Octree *oct)
{
   Photon         p, lastp;
   RREAL          *key;
   OOC_MortonIdx  zIdx, lastzIdx = 0;
   
   rewind(file);
   memset(&lastp, 0, sizeof(lastp));
   
   while (fread(&p, sizeof(p), 1, file) > 0) {
      key = OOC_PhotonKey(&p);
      zIdx = OOC_KEY2MORTON(key, oct);
      
      if (zIdx < lastzIdx) {
         error(INTERNAL, "photons not sorted");
         return -1;
      }
      
      if (zIdx == lastzIdx) {
         sprintf(errmsg, "identical key %021ld "
                 "for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]",
                 zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2], 
                 p.pos [0], p.pos [1], p.pos [2]);
         error(WARNING, errmsg);
      } 
            
      lastzIdx = zIdx;
      memcpy(&lastp, &p, sizeof(p));
   }

   return 0;
}
#endif


void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc)
{
   FILE        *leafFile;
   char        leafFname [1024];
   FVECT       d, octOrg;
   int         i;
   RREAL       octSize = 0;
   
   /* Determine octree size and origin from pmap extent and init octree */
   VCOPY(octOrg, pmap -> minPos);
   VSUB(d, pmap -> maxPos, pmap -> minPos);
   for (i = 0; i < 3; i++)
      if (octSize < d [i])
         octSize = d [i];
         
   if (octSize < FTINY)
      error(INTERNAL, "zero octree size in OOC_BuildPhotonMap");
            
   /* Derive leaf filename from photon map and open file */
   strncpy(leafFname, pmap -> fileName, sizeof(leafFname));
   strncat(leafFname, PMAP_OOC_LEAFSUFFIX, 
           sizeof(leafFname) - strlen(leafFname) - 1);   
   if (!(leafFile = fopen(leafFname, "w+b")))
      error(SYSTEM, "failed to open leaf file in OOC_BuildPhotonMap");

#ifdef DEBUG_OOC
   eputs("Sorting photons by Morton code...\n");
#endif
   
   /* Sort photons in heapfile by Morton code and write to leaf file */
   if (OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE,
                numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey))
      error(INTERNAL, "failed out-of-core photon sort in OOC_BuildPhotonMap");

   /* Init and build octree */   
   OOC_Init(&pmap -> store, sizeof(Photon), octOrg, octSize, OOC_PhotonKey, 
            leafFile);
            
#ifdef DEBUG_OOC
   eputs("Checking leaf file consistency...\n");
   OOC_CheckKeys(leafFile, &pmap -> store);
   
   eputs("Building out-of-core octree...\n");
#endif   
            
   if (!OOC_Build(&pmap -> store, PMAP_OOC_LEAFMAX, PMAP_OOC_MAXDEPTH))
      error(INTERNAL, "failed out-of-core octree build in OOC_BuildPhotonMap");
      
#ifdef DEBUG_OOC
   eputs("Checking out-of-core octree consistency...\n");
   if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), 
                 octOrg, octSize, 0))
      error(INTERNAL, "inconsistent out-of-core octree; Time4Harakiri");
#endif
}


int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out)
{
   return OOC_SaveOctree(&pmap -> store, out);
}


int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *nodeFile)
{
   FILE  *leafFile;
   char  leafFname [1024];
   
   /* Derive leaf filename from photon map and open file */
   strncpy(leafFname, pmap -> fileName, sizeof(leafFname));
   strncat(leafFname, PMAP_OOC_LEAFSUFFIX, 
           sizeof(leafFname) - strlen(leafFname) - 1);   
   
   if (!(leafFile = fopen(leafFname, "r")))
      error(SYSTEM, "failed to open leaf file in OOC_LoadPhotons");

   if (OOC_LoadOctree(&pmap -> store, nodeFile, OOC_PhotonKey, leafFile))
      return -1;
      
   return 0;
}


void OOC_InitFindPhotons (struct PhotonMap *pmap)
{
   if (OOC_InitNearest(&pmap -> squeue, pmap -> maxGather + 1,
                       sizeof(Photon)))
      error(SYSTEM, "can't allocate photon search queue");
}


static void OOC_InitPhotonCache (struct PhotonMap *pmap)
/* Initialise OOC photon cache */
{
   static char warn = 1;
   
   if (!pmap -> store.cache && !pmap -> numDensity) {
      if (pmapCacheSize > 0) {
         const unsigned pageSize = pmapCachePageSize * pmap -> maxGather,
                        numPages = pmapCacheSize / pageSize;
         /* Allocate & init I/O cache in octree */
         pmap -> store.cache = malloc(sizeof(OOC_Cache));
         
         if (!pmap -> store.cache || 
             OOC_CacheInit(pmap -> store.cache, numPages, pageSize, 
                           sizeof(Photon))) {
            error(SYSTEM, "failed OOC photon map cache init");
         }
      }
      else if (warn) {
         error(WARNING, "OOC photon map cache DISABLED");
         warn = 0;
      }
   }   
}


typedef struct {
   const PhotonMap   *pmap;
   const float       *norm;
} OOC_FilterData;


int OOC_FilterPhoton (void *p, void *fd)
/* Filter callback for photon kNN search, used by OOC_FindNearest() */
{
   const Photon         *photon = p;
   const OOC_FilterData *filtData = fd; 
   const PhotonMap      *pmap = filtData -> pmap;

   /* Reject photon if normal faces away (ignored for volume photons) with
    * tolerance to account for perturbation; note photon normal is coded
    * in range [-127,127], hence we factor this in */
   if (filtData -> norm && 
       DOT(filtData->norm, photon->norm) <= PMAP_NORM_TOL * 127 * frandom())
      return 0;
      
   if (isContribPmap(pmap)) {
      /* Lookup in contribution photon map; filter according to emitting
       * light source if contrib list set, else accept all */
       
      if (pmap -> srcContrib) {
         OBJREC *srcMod; 
         const int srcIdx = photonSrcIdx(pmap, photon);
      
         if (srcIdx < 0 || srcIdx >= nsources)
            error(INTERNAL, "invalid light source index in photon map");
      
         srcMod = findmaterial(source [srcIdx].so);

         /* Reject photon if contributions from light source which emitted
          * it are not sought */
         if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
            return 0;
      }

      /* Reject non-caustic photon if lookup for caustic contribs */
      if (pmap -> lookupCaustic && !photon -> caustic)
         return 0;
   }
   
   /* Accept photon */
   return 1;   
}


int OOC_FindPhotons (struct PhotonMap *pmap, const FVECT pos, const FVECT norm)
{
   OOC_SearchFilter  filt;
   OOC_FilterData    filtData;
   float             n [3];
   
   /* Lazily init OOC cache */
   if (!pmap -> store.cache)
      OOC_InitPhotonCache(pmap);
         
   /* Set up filter callback */
   filtData.pmap = pmap;
   if (norm)
      VCOPY(n, norm);
   filtData.norm = norm ? n : NULL;
   filt.data = &filtData;
   filt.func = OOC_FilterPhoton;

   pmap -> maxDist2 = OOC_FindNearest(&pmap -> store, 
                                      OOC_ROOT(&pmap -> store), 0, 
                                      pmap -> store.org, pmap -> store.size,
                                      pos, &filt, &pmap -> squeue, 
                                      pmap -> maxDist2);

   if (pmap -> maxDist2 < 0)
      error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons");
   
   /* Return success or failure (empty queue => none found) */
   return pmap -> squeue.tail ? 0 : -1; 
}


int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, 
                     const FVECT norm, Photon *photon)
{
   OOC_SearchFilter     filt;
   OOC_FilterData       filtData;
   float                n [3], maxDist2;
   
   /* Lazily init OOC cache */
   if (!pmap -> store.cache)
      OOC_InitPhotonCache(pmap);
   
   /* Set up filter callback */
   filtData.pmap = pmap;
   if (norm)
      VCOPY(n, norm);
   filtData.norm = norm ? n : NULL;   
   filt.data = &filtData;
   filt.func = OOC_FilterPhoton;
   
   maxDist2 = OOC_Find1Nearest(&pmap -> store, 
                               OOC_ROOT(&pmap -> store), 0, 
                               pmap -> store.org, pmap -> store.size,
                               pos, &filt, photon, pmap -> maxDist2);
                               
   if (maxDist2 < 0)
      error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon");      
      
   if (maxDist2 >= pmap -> maxDist2)
      /* No photon found => failed */
      return -1;
   else {
      /* Set photon distance => success */
      pmap -> maxDist2 = maxDist2;
      return 0;
   }
}


int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, 
                   Photon *photon)
{
   return OOC_GetData(&pmap -> store, idx, photon);
}


Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
{
   return OOC_GetNearest(squeue, idx);
}


PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap)
{
   return 0;
}
Revision:	1.6
Committed:	Wed Apr 8 15:14:21 2020 UTC (4 years, 1 month ago) by rschregle
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R3
Changes since 1.5:	+24 -12 lines
Log Message:	Fixed est00pid bug in single photon lookups that returned junk when none found, added code to detect and handle (ignore).
#	User	Rev	Content
1	rschregle	1.1	/*
2	rschregle	1.2	======================================================================
3	rschregle	1.1	Photon map interface to out-of-core octree
4
5			Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
6			(c) Lucerne University of Applied Sciences and Arts,
7	rschregle	1.2	supported by the Swiss National Science Foundation (SNSF, #147053)
8			======================================================================
9	rschregle	1.1
10	rschregle	1.6	$Id: pmapooc.c,v 1.5 2018/11/08 00:54:07 greg Exp $
11	rschregle	1.1	*/
12
13
14	rschregle	1.2
15	rschregle	1.1	#include "pmapdata.h" /* Includes pmapooc.h */
16			#include "source.h"
17	greg	1.5	#include "otspecial.h"
18	rschregle	1.1	#include "oocsort.h"
19			#include "oocbuild.h"
20
21
22
23			/* Returns photon position as sorting key for OOC_Octree & friends (notably
24			* for Morton code generation).
25			* !!! Uses type conversion from float via TEMPORARY storage;
26			* !!! THIS IS NOT THREAD SAFE!
27			* !!! RETURNED KEY PERSISTS ONLY IF COPIED BEFORE NEXT CALL! */
28			RREAL OOC_PhotonKey (const void p)
29			{
30			static FVECT photonPos; /* Temp storage for type conversion */
31
32			VCOPY(photonPos, ((Photon*)p) -> pos);
33			return photonPos;
34			}
35
36
37
38			#ifdef DEBUG_OOC
39			static int OOC_CheckKeys (FILE file, const OOC_Octree oct)
40			{
41			Photon p, lastp;
42			RREAL *key;
43			OOC_MortonIdx zIdx, lastzIdx = 0;
44
45			rewind(file);
46			memset(&lastp, 0, sizeof(lastp));
47
48			while (fread(&p, sizeof(p), 1, file) > 0) {
49			key = OOC_PhotonKey(&p);
50			zIdx = OOC_KEY2MORTON(key, oct);
51
52			if (zIdx < lastzIdx) {
53			error(INTERNAL, "photons not sorted");
54			return -1;
55			}
56
57			if (zIdx == lastzIdx) {
58			sprintf(errmsg, "identical key %021ld "
59			"for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]",
60			zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2],
61			p.pos [0], p.pos [1], p.pos [2]);
62			error(WARNING, errmsg);
63			}
64
65			lastzIdx = zIdx;
66			memcpy(&lastp, &p, sizeof(p));
67			}
68
69			return 0;
70			}
71			#endif
72
73
74
75			void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc)
76			{
77			FILE *leafFile;
78			char leafFname [1024];
79			FVECT d, octOrg;
80			int i;
81			RREAL octSize = 0;
82
83			/* Determine octree size and origin from pmap extent and init octree */
84			VCOPY(octOrg, pmap -> minPos);
85			VSUB(d, pmap -> maxPos, pmap -> minPos);
86			for (i = 0; i < 3; i++)
87			if (octSize < d [i])
88			octSize = d [i];
89
90			if (octSize < FTINY)
91			error(INTERNAL, "zero octree size in OOC_BuildPhotonMap");
92
93			/* Derive leaf filename from photon map and open file */
94			strncpy(leafFname, pmap -> fileName, sizeof(leafFname));
95			strncat(leafFname, PMAP_OOC_LEAFSUFFIX,
96			sizeof(leafFname) - strlen(leafFname) - 1);
97			if (!(leafFile = fopen(leafFname, "w+b")))
98			error(SYSTEM, "failed to open leaf file in OOC_BuildPhotonMap");
99
100			#ifdef DEBUG_OOC
101			eputs("Sorting photons by Morton code...\n");
102			#endif
103
104			/* Sort photons in heapfile by Morton code and write to leaf file */
105			if (OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE,
106			numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey))
107			error(INTERNAL, "failed out-of-core photon sort in OOC_BuildPhotonMap");
108
109			/* Init and build octree */
110			OOC_Init(&pmap -> store, sizeof(Photon), octOrg, octSize, OOC_PhotonKey,
111			leafFile);
112
113			#ifdef DEBUG_OOC
114			eputs("Checking leaf file consistency...\n");
115			OOC_CheckKeys(leafFile, &pmap -> store);
116
117			eputs("Building out-of-core octree...\n");
118			#endif
119
120			if (!OOC_Build(&pmap -> store, PMAP_OOC_LEAFMAX, PMAP_OOC_MAXDEPTH))
121			error(INTERNAL, "failed out-of-core octree build in OOC_BuildPhotonMap");
122
123			#ifdef DEBUG_OOC
124			eputs("Checking out-of-core octree consistency...\n");
125			if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store),
126			octOrg, octSize, 0))
127			error(INTERNAL, "inconsistent out-of-core octree; Time4Harakiri");
128			#endif
129			}
130
131
132
133			int OOC_SavePhotons (const struct PhotonMap pmap, FILE out)
134			{
135			return OOC_SaveOctree(&pmap -> store, out);
136			}
137
138
139
140			int OOC_LoadPhotons (struct PhotonMap pmap, FILE nodeFile)
141			{
142			FILE *leafFile;
143			char leafFname [1024];
144
145			/* Derive leaf filename from photon map and open file */
146			strncpy(leafFname, pmap -> fileName, sizeof(leafFname));
147			strncat(leafFname, PMAP_OOC_LEAFSUFFIX,
148			sizeof(leafFname) - strlen(leafFname) - 1);
149
150			if (!(leafFile = fopen(leafFname, "r")))
151			error(SYSTEM, "failed to open leaf file in OOC_LoadPhotons");
152
153			if (OOC_LoadOctree(&pmap -> store, nodeFile, OOC_PhotonKey, leafFile))
154			return -1;
155
156			return 0;
157			}
158
159
160
161			void OOC_InitFindPhotons (struct PhotonMap *pmap)
162			{
163			if (OOC_InitNearest(&pmap -> squeue, pmap -> maxGather + 1,
164			sizeof(Photon)))
165			error(SYSTEM, "can't allocate photon search queue");
166			}
167
168
169
170			static void OOC_InitPhotonCache (struct PhotonMap *pmap)
171			/* Initialise OOC photon cache */
172			{
173			static char warn = 1;
174
175			if (!pmap -> store.cache && !pmap -> numDensity) {
176			if (pmapCacheSize > 0) {
177			const unsigned pageSize = pmapCachePageSize * pmap -> maxGather,
178			numPages = pmapCacheSize / pageSize;
179			/* Allocate & init I/O cache in octree */
180			pmap -> store.cache = malloc(sizeof(OOC_Cache));
181
182			if (!pmap -> store.cache \|\|
183			OOC_CacheInit(pmap -> store.cache, numPages, pageSize,
184			sizeof(Photon))) {
185			error(SYSTEM, "failed OOC photon map cache init");
186			}
187			}
188			else if (warn) {
189			error(WARNING, "OOC photon map cache DISABLED");
190			warn = 0;
191			}
192			}
193			}
194
195
196
197			typedef struct {
198			const PhotonMap *pmap;
199			const float *norm;
200			} OOC_FilterData;
201
202
203
204			int OOC_FilterPhoton (void p, void fd)
205			/* Filter callback for photon kNN search, used by OOC_FindNearest() */
206			{
207			const Photon *photon = p;
208			const OOC_FilterData *filtData = fd;
209			const PhotonMap *pmap = filtData -> pmap;
210
211			/* Reject photon if normal faces away (ignored for volume photons) with
212			* tolerance to account for perturbation; note photon normal is coded
213			* in range [-127,127], hence we factor this in */
214			if (filtData -> norm &&
215			DOT(filtData->norm, photon->norm) <= PMAP_NORM_TOL * 127 * frandom())
216			return 0;
217
218	rschregle	1.3	if (isContribPmap(pmap)) {
219			/* Lookup in contribution photon map; filter according to emitting
220			* light source if contrib list set, else accept all */
221
222			if (pmap -> srcContrib) {
223			OBJREC *srcMod;
224			const int srcIdx = photonSrcIdx(pmap, photon);
225	rschregle	1.1
226	rschregle	1.3	if (srcIdx < 0 \|\| srcIdx >= nsources)
227			error(INTERNAL, "invalid light source index in photon map");
228	rschregle	1.1
229	rschregle	1.3	srcMod = findmaterial(source [srcIdx].so);
230	rschregle	1.1
231	rschregle	1.3	/* Reject photon if contributions from light source which emitted
232			* it are not sought */
233			if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
234			return 0;
235			}
236	rschregle	1.1
237			/* Reject non-caustic photon if lookup for caustic contribs */
238			if (pmap -> lookupCaustic && !photon -> caustic)
239			return 0;
240			}
241
242			/* Accept photon */
243			return 1;
244			}
245
246
247
248	rschregle	1.6	int OOC_FindPhotons (struct PhotonMap *pmap, const FVECT pos, const FVECT norm)
249	rschregle	1.1	{
250			OOC_SearchFilter filt;
251			OOC_FilterData filtData;
252			float n [3];
253
254			/* Lazily init OOC cache */
255			if (!pmap -> store.cache)
256			OOC_InitPhotonCache(pmap);
257
258			/* Set up filter callback */
259			filtData.pmap = pmap;
260	rschregle	1.4	if (norm)
261			VCOPY(n, norm);
262			filtData.norm = norm ? n : NULL;
263	rschregle	1.1	filt.data = &filtData;
264			filt.func = OOC_FilterPhoton;
265
266			pmap -> maxDist2 = OOC_FindNearest(&pmap -> store,
267			OOC_ROOT(&pmap -> store), 0,
268			pmap -> store.org, pmap -> store.size,
269			pos, &filt, &pmap -> squeue,
270			pmap -> maxDist2);
271
272			if (pmap -> maxDist2 < 0)
273			error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons");
274	rschregle	1.6
275			/* Return success or failure (empty queue => none found) */
276			return pmap -> squeue.tail ? 0 : -1;
277	rschregle	1.1	}
278
279
280
281	rschregle	1.6	int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos,
282			const FVECT norm, Photon *photon)
283	rschregle	1.1	{
284			OOC_SearchFilter filt;
285			OOC_FilterData filtData;
286	rschregle	1.6	float n [3], maxDist2;
287	rschregle	1.1
288			/* Lazily init OOC cache */
289			if (!pmap -> store.cache)
290			OOC_InitPhotonCache(pmap);
291
292			/* Set up filter callback */
293			filtData.pmap = pmap;
294	rschregle	1.4	if (norm)
295			VCOPY(n, norm);
296			filtData.norm = norm ? n : NULL;
297	rschregle	1.1	filt.data = &filtData;
298			filt.func = OOC_FilterPhoton;
299
300	rschregle	1.6	maxDist2 = OOC_Find1Nearest(&pmap -> store,
301			OOC_ROOT(&pmap -> store), 0,
302			pmap -> store.org, pmap -> store.size,
303			pos, &filt, photon, pmap -> maxDist2);
304
305			if (maxDist2 < 0)
306			error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon");
307
308			if (maxDist2 >= pmap -> maxDist2)
309			/* No photon found => failed */
310			return -1;
311			else {
312			/* Set photon distance => success */
313			pmap -> maxDist2 = maxDist2;
314			return 0;
315			}
316	rschregle	1.1	}
317
318
319
320			int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx,
321			Photon *photon)
322			{
323			return OOC_GetData(&pmap -> store, idx, photon);
324			}
325
326
327
328			Photon OOC_GetNearestPhoton (const PhotonSearchQueue squeue, PhotonIdx idx)
329			{
330			return OOC_GetNearest(squeue, idx);
331			}
332
333
334
335			PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap)
336			{
337			return 0;
338			}