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.3 2018/01/24 19:39:05 rschregle Exp $
*/


#include "pmapdata.h"   /* Includes pmapooc.h */
#include "source.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;   
}


void 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");
}


void OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, 
                      const FVECT norm, Photon *photon)
{
   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_Find1Nearest(&pmap -> store, 
                                       OOC_ROOT(&pmap -> store), 0, 
                                       pmap -> store.org, pmap -> store.size,
                                       pos, &filt, photon, pmap -> maxDist2);

   if (pmap -> maxDist2 < 0)
      error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon");
}


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