ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapdata.c
Revision: 2.15
Committed: Tue May 17 17:39:47 2016 UTC (8 years ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.14: +383 -504 lines
Log Message:
Initial import of ooC photon map

File Contents

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