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.6 2020/04/08 15:14:21 rschregle 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;

#ifdef DEBUG_OOC
   /* Check octree for consistency */
   if (OOC_Check(
      &pmap -> store, OOC_ROOT(&pmap -> store),
      pmap -> store.org, pmap -> store.size, 0
   ))
      return -1;
#endif

   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.7
Committed:	Mon Mar 22 23:00:00 2021 UTC (3 years, 2 months ago) by rschregle
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, HEAD
Changes since 1.6:	+11 -2 lines
Log Message:	Added optional consistency check of loaded out-of-core photon map in OOC_LoadPhotons() when compiling for debugging mode (DEBUG_OOC).
#	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.7	$Id: pmapooc.c,v 1.6 2020/04/08 15:14:21 rschregle 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	rschregle	1.7
156			#ifdef DEBUG_OOC
157			/* Check octree for consistency */
158			if (OOC_Check(
159			&pmap -> store, OOC_ROOT(&pmap -> store),
160			pmap -> store.org, pmap -> store.size, 0
161			))
162			return -1;
163			#endif
164
165	rschregle	1.1	return 0;
166			}
167
168
169
170			void OOC_InitFindPhotons (struct PhotonMap *pmap)
171			{
172			if (OOC_InitNearest(&pmap -> squeue, pmap -> maxGather + 1,
173			sizeof(Photon)))
174			error(SYSTEM, "can't allocate photon search queue");
175			}
176
177
178
179			static void OOC_InitPhotonCache (struct PhotonMap *pmap)
180			/* Initialise OOC photon cache */
181			{
182			static char warn = 1;
183
184			if (!pmap -> store.cache && !pmap -> numDensity) {
185			if (pmapCacheSize > 0) {
186			const unsigned pageSize = pmapCachePageSize * pmap -> maxGather,
187			numPages = pmapCacheSize / pageSize;
188			/* Allocate & init I/O cache in octree */
189			pmap -> store.cache = malloc(sizeof(OOC_Cache));
190
191			if (!pmap -> store.cache \|\|
192			OOC_CacheInit(pmap -> store.cache, numPages, pageSize,
193			sizeof(Photon))) {
194			error(SYSTEM, "failed OOC photon map cache init");
195			}
196			}
197			else if (warn) {
198			error(WARNING, "OOC photon map cache DISABLED");
199			warn = 0;
200			}
201			}
202			}
203
204
205
206			typedef struct {
207			const PhotonMap *pmap;
208			const float *norm;
209			} OOC_FilterData;
210
211
212
213			int OOC_FilterPhoton (void p, void fd)
214			/* Filter callback for photon kNN search, used by OOC_FindNearest() */
215			{
216			const Photon *photon = p;
217			const OOC_FilterData *filtData = fd;
218			const PhotonMap *pmap = filtData -> pmap;
219
220			/* Reject photon if normal faces away (ignored for volume photons) with
221			* tolerance to account for perturbation; note photon normal is coded
222			* in range [-127,127], hence we factor this in */
223			if (filtData -> norm &&
224			DOT(filtData->norm, photon->norm) <= PMAP_NORM_TOL * 127 * frandom())
225			return 0;
226
227	rschregle	1.3	if (isContribPmap(pmap)) {
228			/* Lookup in contribution photon map; filter according to emitting
229			* light source if contrib list set, else accept all */
230
231			if (pmap -> srcContrib) {
232			OBJREC *srcMod;
233			const int srcIdx = photonSrcIdx(pmap, photon);
234	rschregle	1.1
235	rschregle	1.3	if (srcIdx < 0 \|\| srcIdx >= nsources)
236			error(INTERNAL, "invalid light source index in photon map");
237	rschregle	1.1
238	rschregle	1.3	srcMod = findmaterial(source [srcIdx].so);
239	rschregle	1.1
240	rschregle	1.3	/* Reject photon if contributions from light source which emitted
241			* it are not sought */
242			if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
243			return 0;
244			}
245	rschregle	1.1
246			/* Reject non-caustic photon if lookup for caustic contribs */
247			if (pmap -> lookupCaustic && !photon -> caustic)
248			return 0;
249			}
250
251			/* Accept photon */
252			return 1;
253			}
254
255
256
257	rschregle	1.6	int OOC_FindPhotons (struct PhotonMap *pmap, const FVECT pos, const FVECT norm)
258	rschregle	1.1	{
259			OOC_SearchFilter filt;
260			OOC_FilterData filtData;
261			float n [3];
262
263			/* Lazily init OOC cache */
264			if (!pmap -> store.cache)
265			OOC_InitPhotonCache(pmap);
266
267			/* Set up filter callback */
268			filtData.pmap = pmap;
269	rschregle	1.4	if (norm)
270			VCOPY(n, norm);
271			filtData.norm = norm ? n : NULL;
272	rschregle	1.1	filt.data = &filtData;
273			filt.func = OOC_FilterPhoton;
274
275			pmap -> maxDist2 = OOC_FindNearest(&pmap -> store,
276			OOC_ROOT(&pmap -> store), 0,
277			pmap -> store.org, pmap -> store.size,
278			pos, &filt, &pmap -> squeue,
279			pmap -> maxDist2);
280
281			if (pmap -> maxDist2 < 0)
282			error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons");
283	rschregle	1.6
284			/* Return success or failure (empty queue => none found) */
285			return pmap -> squeue.tail ? 0 : -1;
286	rschregle	1.1	}
287
288
289
290	rschregle	1.6	int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos,
291			const FVECT norm, Photon *photon)
292	rschregle	1.1	{
293			OOC_SearchFilter filt;
294			OOC_FilterData filtData;
295	rschregle	1.6	float n [3], maxDist2;
296	rschregle	1.1
297			/* Lazily init OOC cache */
298			if (!pmap -> store.cache)
299			OOC_InitPhotonCache(pmap);
300
301			/* Set up filter callback */
302			filtData.pmap = pmap;
303	rschregle	1.4	if (norm)
304			VCOPY(n, norm);
305			filtData.norm = norm ? n : NULL;
306	rschregle	1.1	filt.data = &filtData;
307			filt.func = OOC_FilterPhoton;
308
309	rschregle	1.6	maxDist2 = OOC_Find1Nearest(&pmap -> store,
310			OOC_ROOT(&pmap -> store), 0,
311			pmap -> store.org, pmap -> store.size,
312			pos, &filt, photon, pmap -> maxDist2);
313
314			if (maxDist2 < 0)
315			error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon");
316
317			if (maxDist2 >= pmap -> maxDist2)
318			/* No photon found => failed */
319			return -1;
320			else {
321			/* Set photon distance => success */
322			pmap -> maxDist2 = maxDist2;
323			return 0;
324			}
325	rschregle	1.1	}
326
327
328
329			int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx,
330			Photon *photon)
331			{
332			return OOC_GetData(&pmap -> store, idx, photon);
333			}
334
335
336
337			Photon OOC_GetNearestPhoton (const PhotonSearchQueue squeue, PhotonIdx idx)
338			{
339			return OOC_GetNearest(squeue, idx);
340			}
341
342
343
344			PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap)
345			{
346			return 0;
347			}