1 |
#ifndef lint |
2 |
static const char RCSid[] = "$Id: pmutil.c,v 2.5 2020/04/08 15:14:21 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.5 2020/04/08 15:14:21 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 |
/* 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 |
} |
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 |
float r2; |
142 |
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 |
/* 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 |
|
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 |
scalecolor(flux, 2 * (1 - sqn -> dist2 / r2)); |
184 |
#endif |
185 |
addcolor(irrad, flux); |
186 |
} |
187 |
|
188 |
/* Divide by search area PI * r^2, 1 / PI required as ambient |
189 |
normalisation factor */ |
190 |
scalecolor(irrad, 1 / (PI * PI * r2)); |
191 |
|
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 |
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 |
} |
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 |
float r2, gecc2, ph; |
228 |
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 |
/* 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 |
|
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 |
scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2))); |
277 |
return; |
278 |
} |
279 |
#if 0 |
280 |
else |
281 |
/* Apply bias compensation to density estimate */ |
282 |
volumeBiasComp(pmap, ray, irrad); |
283 |
#endif |
284 |
} |