1 |
greg |
2.11 |
#ifndef lint |
2 |
rschregle |
2.23 |
static const char RCSid[] = "$Id: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $"; |
3 |
greg |
2.11 |
#endif |
4 |
rschregle |
2.15 |
|
5 |
greg |
2.1 |
/* |
6 |
rschregle |
2.15 |
========================================================================= |
7 |
|
|
Photon map types and interface to nearest neighbour lookups in underlying |
8 |
|
|
point cloud data structure. |
9 |
|
|
|
10 |
|
|
The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). |
11 |
|
|
This can be overriden with the PMAP_OOC compiletime switch, which enables |
12 |
|
|
an out-of-core octree (see oococt.{h,c}). |
13 |
greg |
2.1 |
|
14 |
|
|
Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) |
15 |
|
|
(c) Fraunhofer Institute for Solar Energy Systems, |
16 |
rschregle |
2.23 |
supported by the German Research Foundation |
17 |
|
|
(DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) |
18 |
rschregle |
2.3 |
(c) Lucerne University of Applied Sciences and Arts, |
19 |
rschregle |
2.23 |
supported by the Swiss National Science Foundation |
20 |
|
|
(SNSF #147053, "Daylight Redirecting Components") |
21 |
rschregle |
2.15 |
========================================================================== |
22 |
greg |
2.1 |
|
23 |
rschregle |
2.23 |
$Id: pmapdata.c,v 2.22 2020/04/08 15:14:21 rschregle Exp $ |
24 |
greg |
2.1 |
*/ |
25 |
|
|
|
26 |
|
|
|
27 |
|
|
|
28 |
rschregle |
2.18 |
#include "pmapdata.h" |
29 |
greg |
2.1 |
#include "pmaprand.h" |
30 |
|
|
#include "pmapmat.h" |
31 |
|
|
#include "otypes.h" |
32 |
|
|
#include "source.h" |
33 |
|
|
#include "rcontrib.h" |
34 |
rschregle |
2.7 |
#include "random.h" |
35 |
greg |
2.1 |
|
36 |
|
|
|
37 |
|
|
|
38 |
|
|
PhotonMap *photonMaps [NUM_PMAP_TYPES] = { |
39 |
|
|
NULL, NULL, NULL, NULL, NULL, NULL |
40 |
|
|
}; |
41 |
|
|
|
42 |
|
|
|
43 |
|
|
|
44 |
rschregle |
2.15 |
/* Include routines to handle underlying point cloud data structure */ |
45 |
|
|
#ifdef PMAP_OOC |
46 |
|
|
#include "pmapooc.c" |
47 |
|
|
#else |
48 |
|
|
#include "pmapkdt.c" |
49 |
|
|
#endif |
50 |
|
|
|
51 |
rschregle |
2.19 |
/* Ambient include/exclude set (from ambient.c) */ |
52 |
|
|
#ifndef MAXASET |
53 |
|
|
#define MAXASET 4095 |
54 |
|
|
#endif |
55 |
|
|
extern OBJECT ambset [MAXASET+1]; |
56 |
|
|
|
57 |
rschregle |
2.21 |
/* Callback to print photon attributes acc. to user defined format */ |
58 |
|
|
int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); |
59 |
|
|
|
60 |
rschregle |
2.15 |
|
61 |
|
|
|
62 |
greg |
2.1 |
void initPhotonMap (PhotonMap *pmap, PhotonMapType t) |
63 |
|
|
/* Init photon map 'n' stuff... */ |
64 |
|
|
{ |
65 |
|
|
if (!pmap) |
66 |
|
|
return; |
67 |
|
|
|
68 |
rschregle |
2.15 |
pmap -> numPhotons = 0; |
69 |
greg |
2.1 |
pmap -> biasCompHist = NULL; |
70 |
|
|
pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; |
71 |
|
|
pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; |
72 |
|
|
pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; |
73 |
|
|
pmap -> gatherTolerance = gatherTolerance; |
74 |
|
|
pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; |
75 |
|
|
pmap -> numDensity = 0; |
76 |
|
|
pmap -> distribRatio = 1; |
77 |
|
|
pmap -> type = t; |
78 |
rschregle |
2.15 |
pmap -> squeue.node = NULL; |
79 |
|
|
pmap -> squeue.len = 0; |
80 |
greg |
2.1 |
|
81 |
|
|
/* Init local RNG state */ |
82 |
|
|
pmap -> randState [0] = 10243; |
83 |
|
|
pmap -> randState [1] = 39829; |
84 |
|
|
pmap -> randState [2] = 9433; |
85 |
|
|
pmapSeed(randSeed, pmap -> randState); |
86 |
|
|
|
87 |
|
|
/* Set up type-specific photon lookup callback */ |
88 |
|
|
pmap -> lookup = pmapLookup [t]; |
89 |
|
|
|
90 |
rschregle |
2.15 |
/* Mark primary photon ray as unused */ |
91 |
|
|
pmap -> lastPrimary.srcIdx = -1; |
92 |
|
|
pmap -> numPrimary = 0; |
93 |
|
|
pmap -> primaries = NULL; |
94 |
|
|
|
95 |
|
|
/* Init storage */ |
96 |
|
|
pmap -> heap = NULL; |
97 |
|
|
pmap -> heapBuf = NULL; |
98 |
|
|
pmap -> heapBufLen = 0; |
99 |
|
|
#ifdef PMAP_OOC |
100 |
|
|
OOC_Null(&pmap -> store); |
101 |
|
|
#else |
102 |
|
|
kdT_Null(&pmap -> store); |
103 |
|
|
#endif |
104 |
greg |
2.1 |
} |
105 |
|
|
|
106 |
|
|
|
107 |
|
|
|
108 |
rschregle |
2.15 |
void initPhotonHeap (PhotonMap *pmap) |
109 |
greg |
2.1 |
{ |
110 |
rschregle |
2.15 |
int fdFlags; |
111 |
greg |
2.1 |
|
112 |
rschregle |
2.15 |
if (!pmap) |
113 |
|
|
error(INTERNAL, "undefined photon map in initPhotonHeap"); |
114 |
greg |
2.1 |
|
115 |
rschregle |
2.15 |
if (!pmap -> heap) { |
116 |
|
|
/* Open heap file */ |
117 |
rschregle |
2.18 |
mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); |
118 |
|
|
if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) |
119 |
rschregle |
2.15 |
error(SYSTEM, "failed opening heap file in initPhotonHeap"); |
120 |
rschregle |
2.18 |
|
121 |
greg |
2.16 |
#ifdef F_SETFL /* XXX is there an alternate needed for Windows? */ |
122 |
rschregle |
2.15 |
fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); |
123 |
|
|
fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); |
124 |
rschregle |
2.18 |
#endif/* ftruncate(fileno(pmap -> heap), 0); */ |
125 |
greg |
2.1 |
} |
126 |
rschregle |
2.15 |
} |
127 |
|
|
|
128 |
|
|
|
129 |
|
|
|
130 |
|
|
void flushPhotonHeap (PhotonMap *pmap) |
131 |
|
|
{ |
132 |
|
|
int fd; |
133 |
|
|
const unsigned long len = pmap -> heapBufLen * sizeof(Photon); |
134 |
|
|
|
135 |
|
|
if (!pmap) |
136 |
|
|
error(INTERNAL, "undefined photon map in flushPhotonHeap"); |
137 |
|
|
|
138 |
rschregle |
2.18 |
if (!pmap -> heap || !pmap -> heapBuf) { |
139 |
|
|
/* Silently ignore undefined heap |
140 |
|
|
error(INTERNAL, "undefined heap in flushPhotonHeap"); */ |
141 |
|
|
return; |
142 |
|
|
} |
143 |
rschregle |
2.15 |
|
144 |
|
|
/* Atomically seek and write block to heap */ |
145 |
|
|
/* !!! Unbuffered I/O via pwrite() avoids potential race conditions |
146 |
|
|
* !!! and buffer corruption which can occur with lseek()/fseek() |
147 |
|
|
* !!! followed by write()/fwrite(). */ |
148 |
|
|
fd = fileno(pmap -> heap); |
149 |
greg |
2.1 |
|
150 |
rschregle |
2.15 |
#ifdef DEBUG_PMAP |
151 |
|
|
sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), |
152 |
|
|
pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon)); |
153 |
|
|
eputs(errmsg); |
154 |
|
|
#endif |
155 |
|
|
|
156 |
|
|
/*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ |
157 |
|
|
if (write(fd, pmap -> heapBuf, len) != len) |
158 |
|
|
error(SYSTEM, "failed append to heap file in flushPhotonHeap"); |
159 |
greg |
2.17 |
|
160 |
rschregle |
2.18 |
#if NIX |
161 |
rschregle |
2.15 |
if (fsync(fd)) |
162 |
|
|
error(SYSTEM, "failed fsync in flushPhotonHeap"); |
163 |
greg |
2.17 |
#endif |
164 |
rschregle |
2.18 |
|
165 |
rschregle |
2.15 |
pmap -> heapBufLen = 0; |
166 |
|
|
} |
167 |
greg |
2.1 |
|
168 |
rschregle |
2.15 |
|
169 |
|
|
|
170 |
rschregle |
2.18 |
#ifdef DEBUG_PMAP |
171 |
rschregle |
2.15 |
static int checkPhotonHeap (FILE *file) |
172 |
|
|
/* Check heap for nonsensical or duplicate photons */ |
173 |
|
|
{ |
174 |
|
|
Photon p, lastp; |
175 |
|
|
int i, dup; |
176 |
|
|
|
177 |
|
|
rewind(file); |
178 |
|
|
memset(&lastp, 0, sizeof(lastp)); |
179 |
greg |
2.1 |
|
180 |
rschregle |
2.15 |
while (fread(&p, sizeof(p), 1, file)) { |
181 |
|
|
dup = 1; |
182 |
|
|
|
183 |
|
|
for (i = 0; i <= 2; i++) { |
184 |
|
|
if (p.pos [i] < thescene.cuorg [i] || |
185 |
|
|
p.pos [i] > thescene.cuorg [i] + thescene.cusize) { |
186 |
|
|
|
187 |
|
|
sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", |
188 |
|
|
p.pos [0], p.pos [1], p.pos [2]); |
189 |
|
|
error(WARNING, errmsg); |
190 |
|
|
} |
191 |
|
|
|
192 |
|
|
dup &= p.pos [i] == lastp.pos [i]; |
193 |
|
|
} |
194 |
|
|
|
195 |
|
|
if (dup) { |
196 |
|
|
sprintf(errmsg, |
197 |
|
|
"consecutive duplicate photon in heap at [%f, %f, %f]\n", |
198 |
|
|
p.pos [0], p.pos [1], p.pos [2]); |
199 |
|
|
error(WARNING, errmsg); |
200 |
|
|
} |
201 |
|
|
} |
202 |
|
|
|
203 |
|
|
return 0; |
204 |
greg |
2.1 |
} |
205 |
rschregle |
2.15 |
#endif |
206 |
greg |
2.1 |
|
207 |
|
|
|
208 |
|
|
|
209 |
rschregle |
2.15 |
int newPhoton (PhotonMap* pmap, const RAY* ray) |
210 |
greg |
2.1 |
{ |
211 |
rschregle |
2.19 |
unsigned i; |
212 |
rschregle |
2.15 |
Photon photon; |
213 |
greg |
2.1 |
COLOR photonFlux; |
214 |
|
|
|
215 |
|
|
/* Account for distribution ratio */ |
216 |
|
|
if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) |
217 |
rschregle |
2.15 |
return -1; |
218 |
greg |
2.1 |
|
219 |
|
|
/* Don't store on sources */ |
220 |
|
|
if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) |
221 |
rschregle |
2.15 |
return -1; |
222 |
greg |
2.1 |
|
223 |
rschregle |
2.20 |
/* Ignore photon if modifier in/outside exclude/include set */ |
224 |
rschregle |
2.19 |
if (ambincl != -1 && ray -> ro && |
225 |
|
|
ambincl != inset(ambset, ray -> ro -> omod)) |
226 |
|
|
return -1; |
227 |
|
|
|
228 |
rschregle |
2.23 |
if (pmapNumROI && pmapROI) { |
229 |
rschregle |
2.19 |
unsigned inROI = 0; |
230 |
rschregle |
2.23 |
FVECT photonDist; |
231 |
rschregle |
2.19 |
|
232 |
rschregle |
2.23 |
/* Store photon if within a region of interest (for ze Ecksperts!) |
233 |
|
|
Note size of spherical ROI is squared. */ |
234 |
|
|
for (i = 0; !inROI && i < pmapNumROI; i++) { |
235 |
|
|
VSUB(photonDist, ray -> rop, pmapROI [i].pos); |
236 |
|
|
|
237 |
|
|
inROI = ( |
238 |
|
|
PMAP_ROI_ISSPHERE(pmapROI + i) |
239 |
|
|
? DOT(photonDist, photonDist) <= pmapROI [i].siz [0] |
240 |
|
|
: fabs(photonDist [0]) <= pmapROI [i].siz [0] && |
241 |
|
|
fabs(photonDist [1]) <= pmapROI [i].siz [1] && |
242 |
|
|
fabs(photonDist [2]) <= pmapROI [i].siz [2] |
243 |
|
|
); |
244 |
|
|
} |
245 |
rschregle |
2.19 |
if (!inROI) |
246 |
|
|
return -1; |
247 |
rschregle |
2.18 |
} |
248 |
rschregle |
2.20 |
|
249 |
rschregle |
2.19 |
/* Adjust flux according to distribution ratio and ray weight */ |
250 |
rschregle |
2.21 |
copycolor(photonFlux, ray -> rcol); |
251 |
|
|
#if 0 |
252 |
|
|
/* Factored out ray -> rweight as deprecated (?) for pmap, and infact |
253 |
|
|
erroneously attenuates volume photon flux based on extinction, |
254 |
|
|
which is already factored in by photonParticipate() */ |
255 |
rschregle |
2.19 |
scalecolor(photonFlux, |
256 |
|
|
ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio |
257 |
|
|
: 1)); |
258 |
rschregle |
2.21 |
#else |
259 |
|
|
scalecolor(photonFlux, |
260 |
|
|
1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1)); |
261 |
|
|
#endif |
262 |
rschregle |
2.19 |
setPhotonFlux(&photon, photonFlux); |
263 |
rschregle |
2.20 |
|
264 |
rschregle |
2.19 |
/* Set photon position and flags */ |
265 |
|
|
VCOPY(photon.pos, ray -> rop); |
266 |
|
|
photon.flags = 0; |
267 |
|
|
photon.caustic = PMAP_CAUSTICRAY(ray); |
268 |
|
|
|
269 |
|
|
/* Set contrib photon's primary ray and subprocess index (the latter |
270 |
|
|
* to linearise the primary ray indices after photon distribution is |
271 |
|
|
* complete). Also set primary ray's source index, thereby marking it |
272 |
|
|
* as used. */ |
273 |
|
|
if (isContribPmap(pmap)) { |
274 |
|
|
photon.primary = pmap -> numPrimary; |
275 |
|
|
photon.proc = PMAP_GETRAYPROC(ray); |
276 |
|
|
pmap -> lastPrimary.srcIdx = ray -> rsrc; |
277 |
|
|
} |
278 |
|
|
else photon.primary = 0; |
279 |
|
|
|
280 |
|
|
/* Set normal */ |
281 |
|
|
for (i = 0; i <= 2; i++) |
282 |
|
|
photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i] |
283 |
|
|
: ray -> ron [i]); |
284 |
rschregle |
2.15 |
|
285 |
rschregle |
2.19 |
if (!pmap -> heapBuf) { |
286 |
|
|
/* Lazily allocate heap buffa */ |
287 |
rschregle |
2.18 |
#if NIX |
288 |
rschregle |
2.19 |
/* Randomise buffa size to temporally decorellate flushes in |
289 |
|
|
* multiprocessing mode */ |
290 |
|
|
srandom(randSeed + getpid()); |
291 |
|
|
pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom()); |
292 |
|
|
#else |
293 |
|
|
/* Randomisation disabled for single processes on WIN; also useful |
294 |
|
|
* for reproducability during debugging */ |
295 |
|
|
pmap -> heapBufSize = PMAP_HEAPBUFSIZE; |
296 |
rschregle |
2.15 |
#endif |
297 |
rschregle |
2.19 |
if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) |
298 |
|
|
error(SYSTEM, "failed heap buffer allocation in newPhoton"); |
299 |
|
|
pmap -> heapBufLen = 0; |
300 |
|
|
} |
301 |
rschregle |
2.15 |
|
302 |
rschregle |
2.19 |
/* Photon initialised; now append to heap buffa */ |
303 |
|
|
memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); |
304 |
|
|
|
305 |
|
|
if (++pmap -> heapBufLen >= pmap -> heapBufSize) |
306 |
|
|
/* Heap buffa full, flush to heap file */ |
307 |
|
|
flushPhotonHeap(pmap); |
308 |
greg |
2.1 |
|
309 |
rschregle |
2.19 |
pmap -> numPhotons++; |
310 |
rschregle |
2.21 |
|
311 |
|
|
/* Print photon attributes */ |
312 |
|
|
if (printPhoton) |
313 |
|
|
/* Non-const kludge */ |
314 |
|
|
printPhoton((RAY*)ray, &photon, pmap); |
315 |
greg |
2.1 |
|
316 |
rschregle |
2.15 |
return 0; |
317 |
greg |
2.1 |
} |
318 |
|
|
|
319 |
|
|
|
320 |
|
|
|
321 |
rschregle |
2.15 |
void buildPhotonMap (PhotonMap *pmap, double *photonFlux, |
322 |
|
|
PhotonPrimaryIdx *primaryOfs, unsigned nproc) |
323 |
|
|
{ |
324 |
|
|
unsigned long n, nCheck = 0; |
325 |
|
|
unsigned i; |
326 |
|
|
Photon *p; |
327 |
|
|
COLOR flux; |
328 |
rschregle |
2.18 |
char nuHeapFname [sizeof(PMAP_TMPFNAME)]; |
329 |
rschregle |
2.15 |
FILE *nuHeap; |
330 |
|
|
/* Need double here to reduce summation errors */ |
331 |
|
|
double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0; |
332 |
|
|
FVECT d; |
333 |
|
|
|
334 |
|
|
if (!pmap) |
335 |
|
|
error(INTERNAL, "undefined photon map in buildPhotonMap"); |
336 |
|
|
|
337 |
|
|
/* Get number of photons from heapfile size */ |
338 |
rschregle |
2.18 |
if (fseek(pmap -> heap, 0, SEEK_END) < 0) |
339 |
|
|
error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap"); |
340 |
rschregle |
2.15 |
pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); |
341 |
|
|
|
342 |
|
|
if (!pmap -> numPhotons) |
343 |
|
|
error(INTERNAL, "empty photon map in buildPhotonMap"); |
344 |
|
|
|
345 |
|
|
if (!pmap -> heap) |
346 |
|
|
error(INTERNAL, "no heap in buildPhotonMap"); |
347 |
|
|
|
348 |
|
|
#ifdef DEBUG_PMAP |
349 |
|
|
eputs("Checking photon heap consistency...\n"); |
350 |
|
|
checkPhotonHeap(pmap -> heap); |
351 |
|
|
|
352 |
|
|
sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); |
353 |
|
|
eputs(errmsg); |
354 |
|
|
#endif |
355 |
rschregle |
2.18 |
|
356 |
rschregle |
2.15 |
/* Allocate heap buffa */ |
357 |
|
|
if (!pmap -> heapBuf) { |
358 |
|
|
pmap -> heapBufSize = PMAP_HEAPBUFSIZE; |
359 |
|
|
pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); |
360 |
|
|
if (!pmap -> heapBuf) |
361 |
|
|
error(SYSTEM, "failed to allocate postprocessed photon heap in" |
362 |
|
|
"buildPhotonMap"); |
363 |
greg |
2.1 |
} |
364 |
|
|
|
365 |
rschregle |
2.15 |
/* We REALLY don't need yet another @%&*! heap just to hold the scaled |
366 |
|
|
* photons, but can't think of a quicker fix... */ |
367 |
rschregle |
2.18 |
mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); |
368 |
|
|
if (!(nuHeap = fopen(nuHeapFname, "w+b"))) |
369 |
rschregle |
2.15 |
error(SYSTEM, "failed to open postprocessed photon heap in " |
370 |
|
|
"buildPhotonMap"); |
371 |
|
|
|
372 |
|
|
rewind(pmap -> heap); |
373 |
|
|
|
374 |
|
|
#ifdef DEBUG_PMAP |
375 |
|
|
eputs("Postprocessing photons...\n"); |
376 |
|
|
#endif |
377 |
greg |
2.1 |
|
378 |
rschregle |
2.18 |
while (!feof(pmap -> heap)) { |
379 |
|
|
#ifdef DEBUG_PMAP |
380 |
|
|
printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap)); |
381 |
|
|
#endif |
382 |
rschregle |
2.15 |
pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), |
383 |
rschregle |
2.18 |
pmap -> heapBufSize, pmap -> heap); |
384 |
|
|
#ifdef DEBUG_PMAP |
385 |
|
|
printf("Got %lu\n", pmap -> heapBufLen); |
386 |
|
|
#endif |
387 |
|
|
|
388 |
|
|
if (ferror(pmap -> heap)) |
389 |
|
|
error(SYSTEM, "failed to read photon heap in buildPhotonMap"); |
390 |
rschregle |
2.15 |
|
391 |
rschregle |
2.18 |
for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { |
392 |
|
|
/* Update min and max pos and set photon flux */ |
393 |
|
|
for (i = 0; i <= 2; i++) { |
394 |
|
|
if (p -> pos [i] < pmap -> minPos [i]) |
395 |
|
|
pmap -> minPos [i] = p -> pos [i]; |
396 |
|
|
else if (p -> pos [i] > pmap -> maxPos [i]) |
397 |
|
|
pmap -> maxPos [i] = p -> pos [i]; |
398 |
|
|
|
399 |
|
|
/* Update centre of gravity with photon position */ |
400 |
|
|
CoG [i] += p -> pos [i]; |
401 |
|
|
} |
402 |
|
|
|
403 |
|
|
if (primaryOfs) |
404 |
|
|
/* Linearise photon primary index from subprocess index using the |
405 |
|
|
* per-subprocess offsets in primaryOfs */ |
406 |
|
|
p -> primary += primaryOfs [p -> proc]; |
407 |
|
|
|
408 |
|
|
/* Scale photon's flux (hitherto normalised to 1 over RGB); in |
409 |
|
|
* case of a contrib photon map, this is done per light source, |
410 |
|
|
* and photonFlux is assumed to be an array */ |
411 |
rschregle |
2.20 |
getPhotonFlux(p, flux); |
412 |
rschregle |
2.18 |
|
413 |
|
|
if (photonFlux) { |
414 |
|
|
scalecolor(flux, photonFlux [isContribPmap(pmap) ? |
415 |
|
|
photonSrcIdx(pmap, p) : 0]); |
416 |
|
|
setPhotonFlux(p, flux); |
417 |
greg |
2.1 |
} |
418 |
rschregle |
2.18 |
|
419 |
|
|
/* Update average photon flux; need a double here */ |
420 |
|
|
addcolor(avgFlux, flux); |
421 |
|
|
} |
422 |
greg |
2.1 |
|
423 |
rschregle |
2.18 |
/* Write modified photons to new heap */ |
424 |
|
|
fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); |
425 |
rschregle |
2.15 |
|
426 |
rschregle |
2.18 |
if (ferror(nuHeap)) |
427 |
|
|
error(SYSTEM, "failed postprocessing photon flux in " |
428 |
|
|
"buildPhotonMap"); |
429 |
rschregle |
2.15 |
|
430 |
|
|
nCheck += pmap -> heapBufLen; |
431 |
|
|
} |
432 |
|
|
|
433 |
|
|
#ifdef DEBUG_PMAP |
434 |
|
|
if (nCheck < pmap -> numPhotons) |
435 |
|
|
error(INTERNAL, "truncated photon heap in buildPhotonMap"); |
436 |
|
|
#endif |
437 |
|
|
|
438 |
|
|
/* Finalise average photon flux */ |
439 |
|
|
scalecolor(avgFlux, 1.0 / pmap -> numPhotons); |
440 |
|
|
copycolor(pmap -> photonFlux, avgFlux); |
441 |
|
|
|
442 |
|
|
/* Average photon positions to get centre of gravity */ |
443 |
|
|
for (i = 0; i < 3; i++) |
444 |
|
|
pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; |
445 |
|
|
|
446 |
|
|
rewind(pmap -> heap); |
447 |
|
|
|
448 |
|
|
/* Compute average photon distance to centre of gravity */ |
449 |
|
|
while (!feof(pmap -> heap)) { |
450 |
|
|
pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), |
451 |
rschregle |
2.18 |
pmap -> heapBufSize, pmap -> heap); |
452 |
rschregle |
2.15 |
|
453 |
rschregle |
2.18 |
for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { |
454 |
|
|
VSUB(d, p -> pos, CoG); |
455 |
|
|
CoGdist += DOT(d, d); |
456 |
|
|
} |
457 |
rschregle |
2.15 |
} |
458 |
|
|
|
459 |
|
|
pmap -> CoGdist = CoGdist /= pmap -> numPhotons; |
460 |
|
|
|
461 |
rschregle |
2.18 |
/* Swap heaps, discarding unscaled photons */ |
462 |
rschregle |
2.15 |
fclose(pmap -> heap); |
463 |
rschregle |
2.18 |
unlink(pmap -> heapFname); |
464 |
rschregle |
2.15 |
pmap -> heap = nuHeap; |
465 |
rschregle |
2.18 |
strcpy(pmap -> heapFname, nuHeapFname); |
466 |
rschregle |
2.15 |
|
467 |
|
|
#ifdef PMAP_OOC |
468 |
|
|
OOC_BuildPhotonMap(pmap, nproc); |
469 |
|
|
#else |
470 |
|
|
kdT_BuildPhotonMap(pmap); |
471 |
|
|
#endif |
472 |
|
|
|
473 |
|
|
/* Trash heap and its buffa */ |
474 |
|
|
free(pmap -> heapBuf); |
475 |
|
|
fclose(pmap -> heap); |
476 |
rschregle |
2.18 |
unlink(pmap -> heapFname); |
477 |
rschregle |
2.15 |
pmap -> heap = NULL; |
478 |
|
|
pmap -> heapBuf = NULL; |
479 |
greg |
2.1 |
} |
480 |
|
|
|
481 |
|
|
|
482 |
|
|
|
483 |
|
|
/* Dynamic max photon search radius increase and reduction factors */ |
484 |
|
|
#define PMAP_MAXDIST_INC 4 |
485 |
|
|
#define PMAP_MAXDIST_DEC 0.9 |
486 |
|
|
|
487 |
|
|
/* Num successful lookups before reducing in max search radius */ |
488 |
|
|
#define PMAP_MAXDIST_CNT 1000 |
489 |
|
|
|
490 |
|
|
/* Threshold below which we assume increasing max radius won't help */ |
491 |
|
|
#define PMAP_SHORT_LOOKUP_THRESH 1 |
492 |
|
|
|
493 |
rschregle |
2.8 |
/* Coefficient for adaptive maximum search radius */ |
494 |
|
|
#define PMAP_MAXDIST_COEFF 100 |
495 |
|
|
|
496 |
greg |
2.1 |
void findPhotons (PhotonMap* pmap, const RAY* ray) |
497 |
|
|
{ |
498 |
|
|
int redo = 0; |
499 |
|
|
|
500 |
rschregle |
2.15 |
if (!pmap -> squeue.len) { |
501 |
greg |
2.1 |
/* Lazy init priority queue */ |
502 |
rschregle |
2.15 |
#ifdef PMAP_OOC |
503 |
|
|
OOC_InitFindPhotons(pmap); |
504 |
|
|
#else |
505 |
|
|
kdT_InitFindPhotons(pmap); |
506 |
|
|
#endif |
507 |
greg |
2.1 |
pmap -> minGathered = pmap -> maxGather; |
508 |
|
|
pmap -> maxGathered = pmap -> minGather; |
509 |
|
|
pmap -> totalGathered = 0; |
510 |
|
|
pmap -> numLookups = pmap -> numShortLookups = 0; |
511 |
|
|
pmap -> shortLookupPct = 0; |
512 |
|
|
pmap -> minError = FHUGE; |
513 |
|
|
pmap -> maxError = -FHUGE; |
514 |
|
|
pmap -> rmsError = 0; |
515 |
rschregle |
2.9 |
/* SQUARED max search radius limit is based on avg photon distance to |
516 |
rschregle |
2.8 |
* centre of gravity, unless fixed by user (maxDistFix > 0) */ |
517 |
rschregle |
2.15 |
pmap -> maxDist0 = pmap -> maxDist2Limit = |
518 |
rschregle |
2.9 |
maxDistFix > 0 ? maxDistFix * maxDistFix |
519 |
rschregle |
2.15 |
: PMAP_MAXDIST_COEFF * pmap -> squeue.len * |
520 |
|
|
pmap -> CoGdist / pmap -> numPhotons; |
521 |
greg |
2.1 |
} |
522 |
|
|
|
523 |
|
|
do { |
524 |
rschregle |
2.15 |
pmap -> squeue.tail = 0; |
525 |
|
|
pmap -> maxDist2 = pmap -> maxDist0; |
526 |
greg |
2.1 |
|
527 |
|
|
/* Search position is ray -> rorg for volume photons, since we have no |
528 |
|
|
intersection point. Normals are ignored -- these are incident |
529 |
rschregle |
2.22 |
directions). */ |
530 |
|
|
/* NOTE: status returned by XXX_FindPhotons() is currently ignored; |
531 |
|
|
if no photons are found, an empty queue is returned under the |
532 |
|
|
assumption all photons are too distant to contribute significant |
533 |
|
|
flux. */ |
534 |
greg |
2.1 |
if (isVolumePmap(pmap)) { |
535 |
rschregle |
2.15 |
#ifdef PMAP_OOC |
536 |
|
|
OOC_FindPhotons(pmap, ray -> rorg, NULL); |
537 |
|
|
#else |
538 |
|
|
kdT_FindPhotons(pmap, ray -> rorg, NULL); |
539 |
|
|
#endif |
540 |
greg |
2.1 |
} |
541 |
|
|
else { |
542 |
rschregle |
2.15 |
#ifdef PMAP_OOC |
543 |
|
|
OOC_FindPhotons(pmap, ray -> rop, ray -> ron); |
544 |
|
|
#else |
545 |
|
|
kdT_FindPhotons(pmap, ray -> rop, ray -> ron); |
546 |
|
|
#endif |
547 |
greg |
2.1 |
} |
548 |
rschregle |
2.3 |
|
549 |
rschregle |
2.15 |
#ifdef PMAP_LOOKUP_INFO |
550 |
|
|
fprintf(stderr, "%d/%d %s photons found within radius %.3f " |
551 |
|
|
"at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, |
552 |
|
|
pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), |
553 |
|
|
ray -> rop [0], ray -> rop [1], ray -> rop [2], |
554 |
|
|
ray -> ro ? ray -> ro -> oname : "<null>"); |
555 |
|
|
#endif |
556 |
|
|
|
557 |
|
|
if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) { |
558 |
greg |
2.1 |
/* Short lookup; too few photons found */ |
559 |
rschregle |
2.15 |
if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { |
560 |
greg |
2.1 |
/* Ignore short lookups which return fewer than |
561 |
|
|
* PMAP_SHORT_LOOKUP_THRESH photons under the assumption there |
562 |
|
|
* really are no photons in the vicinity, and increasing the max |
563 |
|
|
* search radius therefore won't help */ |
564 |
rschregle |
2.8 |
#ifdef PMAP_LOOKUP_WARN |
565 |
greg |
2.1 |
sprintf(errmsg, |
566 |
|
|
"%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", |
567 |
rschregle |
2.15 |
pmap -> squeue.tail, pmap -> squeue.len, |
568 |
|
|
pmapName [pmap -> type], |
569 |
|
|
ray -> rop [0], ray -> rop [1], ray -> rop [2], |
570 |
greg |
2.1 |
ray -> ro ? ray -> ro -> oname : "<null>"); |
571 |
|
|
error(WARNING, errmsg); |
572 |
rschregle |
2.8 |
#endif |
573 |
rschregle |
2.3 |
|
574 |
rschregle |
2.8 |
/* Bail out after warning if maxDist is fixed */ |
575 |
|
|
if (maxDistFix > 0) |
576 |
|
|
return; |
577 |
|
|
|
578 |
rschregle |
2.15 |
if (pmap -> maxDist0 < pmap -> maxDist2Limit) { |
579 |
greg |
2.1 |
/* Increase max search radius if below limit & redo search */ |
580 |
|
|
pmap -> maxDist0 *= PMAP_MAXDIST_INC; |
581 |
rschregle |
2.8 |
#ifdef PMAP_LOOKUP_REDO |
582 |
greg |
2.1 |
redo = 1; |
583 |
rschregle |
2.8 |
#endif |
584 |
|
|
#ifdef PMAP_LOOKUP_WARN |
585 |
greg |
2.1 |
sprintf(errmsg, |
586 |
|
|
redo ? "restarting photon lookup with max radius %.1e" |
587 |
|
|
: "max photon lookup radius adjusted to %.1e", |
588 |
rschregle |
2.15 |
pmap -> maxDist0); |
589 |
greg |
2.1 |
error(WARNING, errmsg); |
590 |
rschregle |
2.8 |
#endif |
591 |
greg |
2.1 |
} |
592 |
rschregle |
2.8 |
#ifdef PMAP_LOOKUP_REDO |
593 |
greg |
2.1 |
else { |
594 |
|
|
sprintf(errmsg, "max photon lookup radius clamped to %.1e", |
595 |
rschregle |
2.15 |
pmap -> maxDist0); |
596 |
greg |
2.1 |
error(WARNING, errmsg); |
597 |
|
|
} |
598 |
rschregle |
2.8 |
#endif |
599 |
greg |
2.1 |
} |
600 |
|
|
|
601 |
|
|
/* Reset successful lookup counter */ |
602 |
|
|
pmap -> numLookups = 0; |
603 |
rschregle |
2.3 |
} |
604 |
greg |
2.1 |
else { |
605 |
rschregle |
2.8 |
/* Bail out after warning if maxDist is fixed */ |
606 |
|
|
if (maxDistFix > 0) |
607 |
|
|
return; |
608 |
|
|
|
609 |
greg |
2.1 |
/* Increment successful lookup counter and reduce max search radius if |
610 |
|
|
* wraparound */ |
611 |
|
|
pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; |
612 |
|
|
if (!pmap -> numLookups) |
613 |
|
|
pmap -> maxDist0 *= PMAP_MAXDIST_DEC; |
614 |
|
|
|
615 |
|
|
redo = 0; |
616 |
|
|
} |
617 |
rschregle |
2.8 |
|
618 |
greg |
2.1 |
} while (redo); |
619 |
|
|
} |
620 |
|
|
|
621 |
|
|
|
622 |
|
|
|
623 |
rschregle |
2.22 |
Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon) |
624 |
rschregle |
2.15 |
{ |
625 |
rschregle |
2.22 |
/* Init (squared) search radius to avg photon dist to centre of gravity */ |
626 |
|
|
float maxDist2_0 = pmap -> CoGdist; |
627 |
|
|
int found = 0; |
628 |
|
|
#ifdef PMAP_LOOKUP_REDO |
629 |
|
|
#define REDO 1 |
630 |
|
|
#else |
631 |
|
|
#define REDO 0 |
632 |
|
|
#endif |
633 |
greg |
2.1 |
|
634 |
rschregle |
2.22 |
do { |
635 |
|
|
pmap -> maxDist2 = maxDist2_0; |
636 |
rschregle |
2.15 |
#ifdef PMAP_OOC |
637 |
rschregle |
2.22 |
found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon); |
638 |
rschregle |
2.15 |
#else |
639 |
rschregle |
2.22 |
found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); |
640 |
|
|
#endif |
641 |
|
|
if (found < 0) { |
642 |
|
|
/* Expand search radius to retry */ |
643 |
|
|
maxDist2_0 *= 2; |
644 |
|
|
#ifdef PMAP_LOOKUP_WARN |
645 |
|
|
sprintf(errmsg, "failed 1-NN photon lookup" |
646 |
|
|
#ifdef PMAP_LOOKUP_REDO |
647 |
|
|
", retrying with search radius %.2f", maxDist2_0 |
648 |
|
|
#endif |
649 |
|
|
); |
650 |
|
|
error(WARNING, errmsg); |
651 |
|
|
#endif |
652 |
|
|
} |
653 |
|
|
} while (REDO && found < 0); |
654 |
|
|
|
655 |
|
|
/* Return photon buffer containing valid photon, else NULL */ |
656 |
|
|
return found < 0 ? NULL : photon; |
657 |
greg |
2.1 |
} |
658 |
|
|
|
659 |
|
|
|
660 |
|
|
|
661 |
rschregle |
2.15 |
void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) |
662 |
greg |
2.1 |
{ |
663 |
rschregle |
2.15 |
#ifdef PMAP_OOC |
664 |
rschregle |
2.18 |
if (OOC_GetPhoton(pmap, idx, photon)) |
665 |
rschregle |
2.15 |
#else |
666 |
|
|
if (kdT_GetPhoton(pmap, idx, photon)) |
667 |
|
|
#endif |
668 |
|
|
error(INTERNAL, "failed photon lookup"); |
669 |
greg |
2.1 |
} |
670 |
|
|
|
671 |
|
|
|
672 |
|
|
|
673 |
rschregle |
2.15 |
Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) |
674 |
|
|
{ |
675 |
|
|
#ifdef PMAP_OOC |
676 |
|
|
return OOC_GetNearestPhoton(squeue, idx); |
677 |
|
|
#else |
678 |
|
|
return kdT_GetNearestPhoton(squeue, idx); |
679 |
|
|
#endif |
680 |
greg |
2.1 |
} |
681 |
|
|
|
682 |
|
|
|
683 |
|
|
|
684 |
rschregle |
2.15 |
PhotonIdx firstPhoton (const PhotonMap *pmap) |
685 |
|
|
{ |
686 |
|
|
#ifdef PMAP_OOC |
687 |
|
|
return OOC_FirstPhoton(pmap); |
688 |
|
|
#else |
689 |
|
|
return kdT_FirstPhoton(pmap); |
690 |
|
|
#endif |
691 |
greg |
2.1 |
} |
692 |
|
|
|
693 |
|
|
|
694 |
|
|
|
695 |
|
|
void deletePhotons (PhotonMap* pmap) |
696 |
|
|
{ |
697 |
rschregle |
2.15 |
#ifdef PMAP_OOC |
698 |
|
|
OOC_Delete(&pmap -> store); |
699 |
|
|
#else |
700 |
|
|
kdT_Delete(&pmap -> store); |
701 |
|
|
#endif |
702 |
|
|
|
703 |
|
|
free(pmap -> squeue.node); |
704 |
greg |
2.1 |
free(pmap -> biasCompHist); |
705 |
|
|
|
706 |
rschregle |
2.15 |
pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = |
707 |
|
|
pmap -> squeue.len = pmap -> squeue.tail = 0; |
708 |
greg |
2.1 |
} |