ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapdata.c
Revision: 2.17
Committed: Thu Sep 29 20:51:07 2016 UTC (7 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.16: +5 -4 lines
Log Message:
Conditional compile of fsync() call for Windows

File Contents

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