ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapdata.c
Revision: 2.16
Committed: Wed Sep 28 22:19:18 2016 UTC (7 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.15: +4 -2 lines
Log Message:
Made fcntl() call compile conditionally for Windows

File Contents

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