1 |
#ifndef lint |
2 |
static const char RCSid[] = "$Id: pmutil.c,v 2.4 2020/01/23 18:27:02 rschregle Exp $"; |
3 |
#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 |
$Id: pmutil.c,v 2.4 2020/01/23 18:27:02 rschregle Exp $ |
16 |
*/ |
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 |
photonDensity, photonDensity |
32 |
}; |
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 |
sprintf(errmsg, "multiple %s photon maps, dropping previous", |
79 |
pmapName [type]); |
80 |
error(WARNING, errmsg); |
81 |
deletePhotons(pmaps [type]); |
82 |
free(pmaps [type]); |
83 |
} |
84 |
pmaps [type] = pm; |
85 |
|
86 |
/* 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 |
|
94 |
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 |
|
106 |
if (pm -> maxGather > pm -> numPhotons) { |
107 |
error(WARNING, "adjusting density estimate bandwidth"); |
108 |
pm -> minGather = pm -> maxGather = pm -> numPhotons; |
109 |
} |
110 |
} |
111 |
} |
112 |
|
113 |
|
114 |
|
115 |
void cleanUpPmaps (PhotonMap **pmaps) |
116 |
{ |
117 |
unsigned t; |
118 |
|
119 |
for (t = 0; t < NUM_PMAP_TYPES; t++) { |
120 |
if (pmaps [t]) { |
121 |
deletePhotons(pmaps [t]); |
122 |
free(pmaps [t]); |
123 |
} |
124 |
} |
125 |
} |
126 |
|
127 |
|
128 |
|
129 |
|
130 |
void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) |
131 |
/* Photon density estimate. Returns irradiance at ray -> rop. */ |
132 |
{ |
133 |
unsigned i; |
134 |
float r2; |
135 |
COLOR flux; |
136 |
Photon *photon; |
137 |
const PhotonSearchQueueNode *sqn; |
138 |
|
139 |
setcolor(irrad, 0, 0, 0); |
140 |
|
141 |
if (!pmap -> maxGather) |
142 |
return; |
143 |
|
144 |
/* Ignore sources */ |
145 |
if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype)) |
146 |
return; |
147 |
|
148 |
findPhotons(pmap, ray); |
149 |
|
150 |
/* Need at least 2 photons */ |
151 |
if (pmap -> squeue.tail < 2) { |
152 |
#ifdef PMAP_NONEFOUND |
153 |
sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", |
154 |
ray -> ro ? ray -> ro -> oname : "<null>", |
155 |
ray -> rop [0], ray -> rop [1], ray -> rop [2]); |
156 |
error(WARNING, errmsg); |
157 |
#endif |
158 |
|
159 |
return; |
160 |
} |
161 |
|
162 |
if (pmap -> minGather == pmap -> maxGather) { |
163 |
/* No bias compensation. Just do a plain vanilla estimate */ |
164 |
sqn = pmap -> squeue.node + 1; |
165 |
|
166 |
/* Average radius^2 between furthest two photons to improve accuracy */ |
167 |
r2 = max(sqn -> dist2, (sqn + 1) -> dist2); |
168 |
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); |
169 |
|
170 |
/* Skip the extra photon */ |
171 |
for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { |
172 |
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); |
173 |
getPhotonFlux(photon, flux); |
174 |
#ifdef PMAP_EPANECHNIKOV |
175 |
/* Apply Epanechnikov kernel to photon flux based on photon dist */ |
176 |
scalecolor(flux, 2 * (1 - sqn -> dist2 / r2)); |
177 |
#endif |
178 |
addcolor(irrad, flux); |
179 |
} |
180 |
|
181 |
/* Divide by search area PI * r^2, 1 / PI required as ambient |
182 |
normalisation factor */ |
183 |
scalecolor(irrad, 1 / (PI * PI * r2)); |
184 |
|
185 |
return; |
186 |
} |
187 |
else |
188 |
/* Apply bias compensation to density estimate */ |
189 |
biasComp(pmap, irrad); |
190 |
} |
191 |
|
192 |
|
193 |
|
194 |
|
195 |
void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad) |
196 |
/* Returns precomputed photon density estimate at ray -> rop. */ |
197 |
{ |
198 |
Photon p; |
199 |
|
200 |
setcolor(irrad, 0, 0, 0); |
201 |
|
202 |
/* Ignore sources */ |
203 |
if (r -> ro && islight(objptr(r -> ro -> omod) -> otype)) |
204 |
return; |
205 |
|
206 |
if (find1Photon(preCompPmap, r, &p)) |
207 |
/* p contains a found photon, so get its irradiance, otherwise it |
208 |
* remains zero under the assumption all photons are too distant |
209 |
* to contribute significantly */ |
210 |
getPhotonFlux(&p, irrad); |
211 |
} |
212 |
|
213 |
|
214 |
|
215 |
|
216 |
void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) |
217 |
/* Photon volume density estimate. Returns irradiance at ray -> rop. */ |
218 |
{ |
219 |
unsigned i; |
220 |
float r2, gecc2, ph; |
221 |
COLOR flux; |
222 |
Photon *photon; |
223 |
const PhotonSearchQueueNode *sqn; |
224 |
|
225 |
setcolor(irrad, 0, 0, 0); |
226 |
|
227 |
if (!pmap -> maxGather) |
228 |
return; |
229 |
|
230 |
findPhotons(pmap, ray); |
231 |
|
232 |
/* Need at least 2 photons */ |
233 |
if (pmap -> squeue.tail < 2) |
234 |
return; |
235 |
|
236 |
#if 0 |
237 |
/* Volume biascomp disabled (probably redundant) */ |
238 |
if (pmap -> minGather == pmap -> maxGather) |
239 |
#endif |
240 |
{ |
241 |
/* No bias compensation. Just do a plain vanilla estimate */ |
242 |
gecc2 = ray -> gecc * ray -> gecc; |
243 |
sqn = pmap -> squeue.node + 1; |
244 |
|
245 |
/* Average radius^2 between furthest two photons to improve accuracy */ |
246 |
r2 = max(sqn -> dist2, (sqn + 1) -> dist2); |
247 |
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); |
248 |
|
249 |
/* Skip the extra photon */ |
250 |
for (i = 1; i < pmap -> squeue.tail; i++, sqn++) { |
251 |
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); |
252 |
|
253 |
/* Compute phase function for inscattering from photon */ |
254 |
if (gecc2 <= FTINY) |
255 |
ph = 1; |
256 |
else { |
257 |
ph = DOT(ray -> rdir, photon -> norm) / 127; |
258 |
ph = 1 + gecc2 - 2 * ray -> gecc * ph; |
259 |
ph = (1 - gecc2) / (ph * sqrt(ph)); |
260 |
} |
261 |
|
262 |
getPhotonFlux(photon, flux); |
263 |
scalecolor(flux, ph); |
264 |
addcolor(irrad, flux); |
265 |
} |
266 |
|
267 |
/* Divide by search volume 4 / 3 * PI * r^3 and phase function |
268 |
normalization factor 1 / (4 * PI) */ |
269 |
scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2))); |
270 |
return; |
271 |
} |
272 |
#if 0 |
273 |
else |
274 |
/* Apply bias compensation to density estimate */ |
275 |
volumeBiasComp(pmap, ray, irrad); |
276 |
#endif |
277 |
} |