src/rt/pmutil.c

#ifndef lint
static const char RCSid[] = "$Id: pmutil.c,v 2.2 2018/01/24 19:39:05 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.2 2018/01/24 19:39:05 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]) {
            deletePhotons(pmaps [type]);
            free(pmaps [type]);
         }
         pmaps [type] = pm;
         
         /* Check for invalid density estimate bandwidth */                            
         if (pm -> maxGather > pm -> numPhotons) {
            error(WARNING, "adjusting density estimate bandwidth");
            pm -> minGather = pm -> maxGather = pm -> numPhotons;
         }
      }
}


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;
      
   find1Photon(preCompPmap, r, &p);
   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.3
Committed:	Fri Feb 9 14:57:42 2018 UTC (7 years, 8 months ago) by rschregle
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2
Changes since 2.2:	+13 -13 lines
Log Message:	r -> r2 in photonDensity() to emphasise squared radius
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	rschregle	2.3	static const char RCSid[] = "$Id: pmutil.c,v 2.2 2018/01/24 19:39:05 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.3	$Id: pmutil.c,v 2.2 2018/01/24 19:39:05 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
60			for (t = 0; t < NUM_PMAP_TYPES; t++)
61			if (setPmapParam(&pm, parm + t)) {
62			/* 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
71			/* 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			deletePhotons(pmaps [type]);
79			free(pmaps [type]);
80			}
81			pmaps [type] = pm;
82
83			/* Check for invalid density estimate bandwidth */
84			if (pm -> maxGather > pm -> numPhotons) {
85			error(WARNING, "adjusting density estimate bandwidth");
86			pm -> minGather = pm -> maxGather = pm -> numPhotons;
87			}
88			}
89			}
90
91
92
93			void cleanUpPmaps (PhotonMap **pmaps)
94			{
95			unsigned t;
96
97			for (t = 0; t < NUM_PMAP_TYPES; t++) {
98			if (pmaps [t]) {
99			deletePhotons(pmaps [t]);
100			free(pmaps [t]);
101			}
102			}
103			}
104
105
106
107
108			void photonDensity (PhotonMap pmap, RAY ray, COLOR irrad)
109			/* Photon density estimate. Returns irradiance at ray -> rop. */
110			{
111			unsigned i;
112	rschregle	2.3	float r2;
113	greg	2.1	COLOR flux;
114			Photon *photon;
115			const PhotonSearchQueueNode *sqn;
116
117			setcolor(irrad, 0, 0, 0);
118
119			if (!pmap -> maxGather)
120			return;
121
122			/* Ignore sources */
123			if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype))
124			return;
125
126			findPhotons(pmap, ray);
127
128			/* Need at least 2 photons */
129			if (pmap -> squeue.tail < 2) {
130			#ifdef PMAP_NONEFOUND
131			sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
132			ray -> ro ? ray -> ro -> oname : "<null>",
133			ray -> rop [0], ray -> rop [1], ray -> rop [2]);
134			error(WARNING, errmsg);
135			#endif
136
137			return;
138			}
139
140			if (pmap -> minGather == pmap -> maxGather) {
141			/* No bias compensation. Just do a plain vanilla estimate */
142			sqn = pmap -> squeue.node + 1;
143
144	rschregle	2.3	/* Average radius^2 between furthest two photons to improve accuracy */
145			r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
146			r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
147	greg	2.1
148			/* Skip the extra photon */
149			for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
150			photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
151			getPhotonFlux(photon, flux);
152			#ifdef PMAP_EPANECHNIKOV
153			/* Apply Epanechnikov kernel to photon flux based on photon dist */
154	rschregle	2.3	scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
155	greg	2.1	#endif
156			addcolor(irrad, flux);
157			}
158
159			/* Divide by search area PI * r^2, 1 / PI required as ambient
160			normalisation factor */
161	rschregle	2.3	scalecolor(irrad, 1 / (PI * PI * r2));
162	greg	2.1
163			return;
164			}
165			else
166			/* Apply bias compensation to density estimate */
167			biasComp(pmap, irrad);
168			}
169
170
171
172
173			void photonPreCompDensity (PhotonMap pmap, RAY r, COLOR irrad)
174			/* Returns precomputed photon density estimate at ray -> rop. */
175			{
176			Photon p;
177
178			setcolor(irrad, 0, 0, 0);
179
180			/* Ignore sources */
181			if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
182			return;
183
184			find1Photon(preCompPmap, r, &p);
185			getPhotonFlux(&p, irrad);
186			}
187
188
189
190
191			void volumePhotonDensity (PhotonMap pmap, RAY ray, COLOR irrad)
192			/* Photon volume density estimate. Returns irradiance at ray -> rop. */
193			{
194			unsigned i;
195	rschregle	2.3	float r2, gecc2, ph;
196	greg	2.1	COLOR flux;
197			Photon *photon;
198			const PhotonSearchQueueNode *sqn;
199
200			setcolor(irrad, 0, 0, 0);
201
202			if (!pmap -> maxGather)
203			return;
204
205			findPhotons(pmap, ray);
206
207			/* Need at least 2 photons */
208			if (pmap -> squeue.tail < 2)
209			return;
210
211			#if 0
212			/* Volume biascomp disabled (probably redundant) */
213			if (pmap -> minGather == pmap -> maxGather)
214			#endif
215			{
216			/* No bias compensation. Just do a plain vanilla estimate */
217			gecc2 = ray -> gecc * ray -> gecc;
218			sqn = pmap -> squeue.node + 1;
219
220	rschregle	2.3	/* Average radius^2 between furthest two photons to improve accuracy */
221			r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
222			r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
223	greg	2.1
224			/* Skip the extra photon */
225			for (i = 1; i < pmap -> squeue.tail; i++, sqn++) {
226			photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
227
228			/* Compute phase function for inscattering from photon */
229			if (gecc2 <= FTINY)
230			ph = 1;
231			else {
232			ph = DOT(ray -> rdir, photon -> norm) / 127;
233			ph = 1 + gecc2 - 2 * ray -> gecc * ph;
234			ph = (1 - gecc2) / (ph * sqrt(ph));
235			}
236
237			getPhotonFlux(photon, flux);
238			scalecolor(flux, ph);
239			addcolor(irrad, flux);
240			}
241
242			/* Divide by search volume 4 / 3 * PI * r^3 and phase function
243			normalization factor 1 / (4 * PI) */
244	rschregle	2.3	scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
245	greg	2.1	return;
246			}
247			#if 0
248			else
249			/* Apply bias compensation to density estimate */
250			volumeBiasComp(pmap, ray, irrad);
251			#endif
252			}