src/rt/pmutil.c

#ifndef lint
static const char RCSid[] = "$Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $";
#endif

/* 
   ======================================================================
   Photon map utilities

   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: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $
*/

#include "pmap.h"
#include "pmapio.h"
#include "pmapbias.h"
#include "otypes.h"
#include <sys/stat.h>


extern char *octname;


/* Photon map lookup functions per type */
void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
   photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
   photonDensity, photonDensity
};


void colorNorm (COLOR c)
/* Normalise colour channels to average of 1 */
{
   const float avg = colorAvg(c);
   
   if (!avg) 
      return;
      
   c [0] /= avg;
   c [1] /= avg;
   c [2] /= avg;
}


void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm)
{
   unsigned t;
   struct stat octstat, pmstat;
   PhotonMap *pm;
   PhotonMapType type;

   for (t = 0; t < NUM_PMAP_TYPES; t++)
      if (setPmapParam(&pm, parm + t)) {
         /* Check if photon map newer than octree */
         if (pm -> fileName && octname &&
             !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) && 
             octstat.st_mtime > pmstat.st_mtime) {
            sprintf(errmsg, "photon map in file %s may be stale", 
                    pm -> fileName);
            error(USER, errmsg);
         }

         /* Load photon map from file and get its type */
         if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE)
            error(USER, "failed loading photon map");
            
         /* Assign to appropriate photon map type (deleting previously
          * loaded photon map of same type if necessary) */
         if (pmaps [type]) {
            sprintf(errmsg, "multiple %s photon maps, dropping previous",
                    pmapName [type]);
            error(WARNING, errmsg);
            deletePhotons(pmaps [type]);
            free(pmaps [type]);
         }
         pmaps [type] = pm;

         /* Check for valid density estimate bandwidths */
         if ((pm -> minGather > 1 || pm -> maxGather > 1) &&
             (type == PMAP_TYPE_PRECOMP)) {
            /* Force bwidth to 1 for precomputed pmap */
            error(WARNING, "ignoring bandwidth for precomp photon map");
            pm -> minGather = pm -> maxGather = 1;
         }

         if ((pm -> maxGather > pm -> minGather) && 
             (type == PMAP_TYPE_VOLUME)) {            
            /* Biascomp for volume pmaps (see volumePhotonDensity() below) 
               is considered redundant, and there's probably no point in 
               recovering by using the lower bandwidth, since it's probably
               not what the user wants, so bail out. */
            sprintf(errmsg, 
                    "bias compensation is not available with %s photon maps",
                    pmapName [type]);
            error(USER, errmsg);
         }

         if (pm -> maxGather > pm -> numPhotons) {
            /* Clamp lookup bandwidth to total number of photons (minus one,
               since density estimate gets extra photon to obtain averaged 
               radius) */
            sprintf(
               errmsg, "clamping density estimate bandwidth to %ld",
               pm -> numPhotons
            );
            error(WARNING, errmsg);
            pm -> minGather = pm -> maxGather = pm -> numPhotons - 1;
         }
      }
}


void cleanUpPmaps (PhotonMap **pmaps)
{
   unsigned t;
   
   for (t = 0; t < NUM_PMAP_TYPES; t++) {
      if (pmaps [t]) {
         deletePhotons(pmaps [t]);
         free(pmaps [t]);
      }
   }
}


void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Photon density estimate. Returns irradiance at ray -> rop. */
{
   unsigned                      i;
   float                         r2;
   COLOR                         flux;
   Photon                        *photon;
   const PhotonSearchQueueNode   *sqn;
 
   setcolor(irrad, 0, 0, 0);

   if (!pmap -> maxGather) 
      return;
      
   /* Ignore sources */
   if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype)) 
      return;
         
   findPhotons(pmap, ray);
   
   /* Need at least 2 photons */
   if (pmap -> squeue.tail < 2) {
#ifdef PMAP_NONEFOUND   
      sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", 
              ray -> ro ? ray -> ro -> oname : "<null>",
              ray -> rop [0], ray -> rop [1], ray -> rop [2]);
      error(WARNING, errmsg);
#endif      

      return;
   }

   if (pmap -> minGather == pmap -> maxGather) {
      /* No bias compensation. Just do a plain vanilla estimate */
      sqn = pmap -> squeue.node + 1;
      
      /* Average radius^2 between furthest two photons to improve accuracy */
      r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
      r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
      
      /* Skip the extra photon */
      for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
         photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
         getPhotonFlux(photon, flux);         
#ifdef PMAP_EPANECHNIKOV
         /* Apply Epanechnikov kernel to photon flux based on photon dist */
         scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
#endif   
         addcolor(irrad, flux);
      }
      
      /* Divide by search area PI * r^2, 1 / PI required as ambient 
         normalisation factor */         
      scalecolor(irrad, 1 / (PI * PI * r2)); 
      
      return;
   }
   else 
      /* Apply bias compensation to density estimate */
      biasComp(pmap, irrad);
}


void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
/* Returns precomputed photon density estimate at ray -> rop. */
{
   Photon p;
   
   setcolor(irrad, 0, 0, 0);

   /* Ignore sources */
   if (r -> ro && islight(objptr(r -> ro -> omod) -> otype)) 
      return;
      
   if (find1Photon(preCompPmap, r, &p))
      /* p contains a found photon, so get its irradiance, otherwise it
       * remains zero under the assumption all photons are too distant
       * to contribute significantly */ 
      getPhotonFlux(&p, irrad);
}


void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Photon volume density estimate. Returns irradiance at ray -> rop. */
{
   unsigned                      i;
   float                         r2, gecc2, ph;
   COLOR                         flux;
   Photon                        *photon;
   const PhotonSearchQueueNode   *sqn;

   setcolor(irrad, 0, 0, 0);
   
   if (!pmap -> maxGather) 
      return;
      
   findPhotons(pmap, ray);
   
   /* Need at least 2 photons */
   if (pmap -> squeue.tail < 2) 
      return;

#if 0      
   /* Volume biascomp disabled (probably redundant) */
   if (pmap -> minGather == pmap -> maxGather)
#endif   
   {
      /* No bias compensation. Just do a plain vanilla estimate */
      gecc2 = ray -> gecc * ray -> gecc;
      sqn = pmap -> squeue.node + 1;
      
      /* Average radius^2 between furthest two photons to improve accuracy */      
      r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
      r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
      
      /* Skip the extra photon */
      for (i = 1; i < pmap -> squeue.tail; i++, sqn++) {
         photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
         
         /* Compute phase function for inscattering from photon */
         if (gecc2 <= FTINY) 
            ph = 1;
         else {
            ph = DOT(ray -> rdir, photon -> norm) / 127;
            ph = 1 + gecc2 - 2 * ray -> gecc * ph;
            ph = (1 - gecc2) / (ph * sqrt(ph));
         }
         
         getPhotonFlux(photon, flux);
         scalecolor(flux, ph);
         addcolor(irrad, flux);
      }
      
      /* Divide by search volume 4 / 3 * PI * r^3 and phase function
         normalization factor 1 / (4 * PI) */
      scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
      return;
   }
#if 0
   else 
      /* Apply bias compensation to density estimate */
      volumeBiasComp(pmap, ray, irrad);
#endif      
}
Revision:	2.6
Committed:	Tue Mar 23 00:07:13 2021 UTC (3 years, 1 month ago) by rschregle
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, HEAD
Changes since 2.5:	+18 -11 lines
Log Message:	Fixed off-by-one bug in photon density estimate bandwidth clamping if exceeds number of photons.
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	rschregle	2.6	static const char RCSid[] = "$Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $";
3	greg	2.1	#endif
4
5			/*
6			======================================================================
7			Photon map utilities
8
9			Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
10			(c) Fraunhofer Institute for Solar Energy Systems,
11			(c) Lucerne University of Applied Sciences and Arts,
12			supported by the Swiss National Science Foundation (SNSF, #147053)
13			======================================================================
14
15	rschregle	2.6	$Id: pmutil.c,v 2.5 2020/04/08 15:14:21 rschregle Exp $
16	greg	2.1	*/
17
18			#include "pmap.h"
19			#include "pmapio.h"
20			#include "pmapbias.h"
21			#include "otypes.h"
22			#include <sys/stat.h>
23
24
25			extern char *octname;
26
27
28			/* Photon map lookup functions per type */
29			void (pmapLookup [NUM_PMAP_TYPES])(PhotonMap, RAY*, COLOR) = {
30			photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
31	rschregle	2.2	photonDensity, photonDensity
32	greg	2.1	};
33
34
35
36
37			void colorNorm (COLOR c)
38			/* Normalise colour channels to average of 1 */
39			{
40			const float avg = colorAvg(c);
41
42			if (!avg)
43			return;
44
45			c [0] /= avg;
46			c [1] /= avg;
47			c [2] /= avg;
48			}
49
50
51
52
53			void loadPmaps (PhotonMap *pmaps, const PhotonMapParams parm)
54			{
55			unsigned t;
56			struct stat octstat, pmstat;
57			PhotonMap *pm;
58			PhotonMapType type;
59	rschregle	2.6
60	greg	2.1	for (t = 0; t < NUM_PMAP_TYPES; t++)
61	rschregle	2.6	if (setPmapParam(&pm, parm + t)) {
62	greg	2.1	/* Check if photon map newer than octree */
63			if (pm -> fileName && octname &&
64			!stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) &&
65			octstat.st_mtime > pmstat.st_mtime) {
66			sprintf(errmsg, "photon map in file %s may be stale",
67			pm -> fileName);
68			error(USER, errmsg);
69			}
70	rschregle	2.6
71	greg	2.1	/* Load photon map from file and get its type */
72			if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE)
73			error(USER, "failed loading photon map");
74
75			/* Assign to appropriate photon map type (deleting previously
76			* loaded photon map of same type if necessary) */
77			if (pmaps [type]) {
78	rschregle	2.4	sprintf(errmsg, "multiple %s photon maps, dropping previous",
79			pmapName [type]);
80			error(WARNING, errmsg);
81	greg	2.1	deletePhotons(pmaps [type]);
82			free(pmaps [type]);
83			}
84			pmaps [type] = pm;
85	rschregle	2.6
86	rschregle	2.4	/* Check for valid density estimate bandwidths */
87			if ((pm -> minGather > 1 \|\| pm -> maxGather > 1) &&
88			(type == PMAP_TYPE_PRECOMP)) {
89			/* Force bwidth to 1 for precomputed pmap */
90			error(WARNING, "ignoring bandwidth for precomp photon map");
91			pm -> minGather = pm -> maxGather = 1;
92			}
93	rschregle	2.6
94	rschregle	2.4	if ((pm -> maxGather > pm -> minGather) &&
95			(type == PMAP_TYPE_VOLUME)) {
96			/* Biascomp for volume pmaps (see volumePhotonDensity() below)
97			is considered redundant, and there's probably no point in
98			recovering by using the lower bandwidth, since it's probably
99			not what the user wants, so bail out. */
100			sprintf(errmsg,
101			"bias compensation is not available with %s photon maps",
102			pmapName [type]);
103			error(USER, errmsg);
104			}
105	rschregle	2.6
106	greg	2.1	if (pm -> maxGather > pm -> numPhotons) {
107	rschregle	2.6	/* Clamp lookup bandwidth to total number of photons (minus one,
108			since density estimate gets extra photon to obtain averaged
109			radius) */
110			sprintf(
111			errmsg, "clamping density estimate bandwidth to %ld",
112			pm -> numPhotons
113			);
114			error(WARNING, errmsg);
115			pm -> minGather = pm -> maxGather = pm -> numPhotons - 1;
116			}
117	greg	2.1	}
118			}
119
120
121
122			void cleanUpPmaps (PhotonMap **pmaps)
123			{
124			unsigned t;
125
126			for (t = 0; t < NUM_PMAP_TYPES; t++) {
127			if (pmaps [t]) {
128			deletePhotons(pmaps [t]);
129			free(pmaps [t]);
130			}
131			}
132			}
133
134
135
136
137			void photonDensity (PhotonMap pmap, RAY ray, COLOR irrad)
138			/* Photon density estimate. Returns irradiance at ray -> rop. */
139			{
140			unsigned i;
141	rschregle	2.3	float r2;
142	greg	2.1	COLOR flux;
143			Photon *photon;
144			const PhotonSearchQueueNode *sqn;
145
146			setcolor(irrad, 0, 0, 0);
147
148			if (!pmap -> maxGather)
149			return;
150
151			/* Ignore sources */
152			if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype))
153			return;
154
155			findPhotons(pmap, ray);
156
157			/* Need at least 2 photons */
158			if (pmap -> squeue.tail < 2) {
159			#ifdef PMAP_NONEFOUND
160			sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
161			ray -> ro ? ray -> ro -> oname : "<null>",
162			ray -> rop [0], ray -> rop [1], ray -> rop [2]);
163			error(WARNING, errmsg);
164			#endif
165
166			return;
167			}
168
169			if (pmap -> minGather == pmap -> maxGather) {
170			/* No bias compensation. Just do a plain vanilla estimate */
171			sqn = pmap -> squeue.node + 1;
172
173	rschregle	2.3	/* Average radius^2 between furthest two photons to improve accuracy */
174			r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
175			r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
176	greg	2.1
177			/* Skip the extra photon */
178			for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
179			photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
180			getPhotonFlux(photon, flux);
181			#ifdef PMAP_EPANECHNIKOV
182			/* Apply Epanechnikov kernel to photon flux based on photon dist */
183	rschregle	2.3	scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
184	greg	2.1	#endif
185			addcolor(irrad, flux);
186			}
187
188			/* Divide by search area PI * r^2, 1 / PI required as ambient
189			normalisation factor */
190	rschregle	2.3	scalecolor(irrad, 1 / (PI * PI * r2));
191	greg	2.1
192			return;
193			}
194			else
195			/* Apply bias compensation to density estimate */
196			biasComp(pmap, irrad);
197			}
198
199
200
201
202			void photonPreCompDensity (PhotonMap pmap, RAY r, COLOR irrad)
203			/* Returns precomputed photon density estimate at ray -> rop. */
204			{
205			Photon p;
206
207			setcolor(irrad, 0, 0, 0);
208
209			/* Ignore sources */
210			if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
211			return;
212
213	rschregle	2.5	if (find1Photon(preCompPmap, r, &p))
214			/* p contains a found photon, so get its irradiance, otherwise it
215			* remains zero under the assumption all photons are too distant
216			* to contribute significantly */
217			getPhotonFlux(&p, irrad);
218	greg	2.1	}
219
220
221
222
223			void volumePhotonDensity (PhotonMap pmap, RAY ray, COLOR irrad)
224			/* Photon volume density estimate. Returns irradiance at ray -> rop. */
225			{
226			unsigned i;
227	rschregle	2.3	float r2, gecc2, ph;
228	greg	2.1	COLOR flux;
229			Photon *photon;
230			const PhotonSearchQueueNode *sqn;
231
232			setcolor(irrad, 0, 0, 0);
233
234			if (!pmap -> maxGather)
235			return;
236
237			findPhotons(pmap, ray);
238
239			/* Need at least 2 photons */
240			if (pmap -> squeue.tail < 2)
241			return;
242
243			#if 0
244			/* Volume biascomp disabled (probably redundant) */
245			if (pmap -> minGather == pmap -> maxGather)
246			#endif
247			{
248			/* No bias compensation. Just do a plain vanilla estimate */
249			gecc2 = ray -> gecc * ray -> gecc;
250			sqn = pmap -> squeue.node + 1;
251
252	rschregle	2.3	/* Average radius^2 between furthest two photons to improve accuracy */
253			r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
254			r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
255	greg	2.1
256			/* Skip the extra photon */
257			for (i = 1; i < pmap -> squeue.tail; i++, sqn++) {
258			photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
259
260			/* Compute phase function for inscattering from photon */
261			if (gecc2 <= FTINY)
262			ph = 1;
263			else {
264			ph = DOT(ray -> rdir, photon -> norm) / 127;
265			ph = 1 + gecc2 - 2 * ray -> gecc * ph;
266			ph = (1 - gecc2) / (ph * sqrt(ph));
267			}
268
269			getPhotonFlux(photon, flux);
270			scalecolor(flux, ph);
271			addcolor(irrad, flux);
272			}
273
274			/* Divide by search volume 4 / 3 * PI * r^3 and phase function
275			normalization factor 1 / (4 * PI) */
276	rschregle	2.3	scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
277	greg	2.1	return;
278			}
279			#if 0
280			else
281			/* Apply bias compensation to density estimate */
282			volumeBiasComp(pmap, ray, irrad);
283			#endif
284			}