ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapcontrib.c
Revision: 2.13
Committed: Thu Sep 29 21:51:58 2016 UTC (8 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.12: +3 -206 lines
Log Message:
Separated creation from use of contribution photon maps

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: pmapcontrib.c,v 2.12 2016/05/17 17:39:47 rschregle Exp $";
3 #endif
4
5 /*
6 ======================================================================
7 Photon map support for building light source contributions
8
9 Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
10 (c) Lucerne University of Applied Sciences and Arts,
11 supported by the Swiss National Science Foundation (SNSF, #147053)
12 ======================================================================
13
14 $Id: pmapcontrib.c,v 2.12 2016/05/17 17:39:47 rschregle Exp $
15 */
16
17
18 #include "pmapcontrib.h"
19 #include "pmapmat.h"
20 #include "pmapsrc.h"
21 #include "pmaprand.h"
22 #include "pmapio.h"
23 #include "pmapdiag.h"
24 #include "rcontrib.h"
25 #include "otypes.h"
26 #include <sys/mman.h>
27 #include <sys/wait.h>
28
29
30
31 static PhotonPrimaryIdx newPhotonPrimary (PhotonMap *pmap,
32 const RAY *primRay,
33 FILE *primHeap)
34 /* Add primary ray for emitted photon and save light source index, origin on
35 * source, and emitted direction; used by contrib photons. The current
36 * primary is stored in pmap -> lastPrimary. If the previous primary
37 * contributed photons (has srcIdx >= 0), it's appended to primHeap. If
38 * primRay == NULL, the current primary is still flushed, but no new primary
39 * is set. Returns updated primary counter pmap -> numPrimary. */
40 {
41 if (!pmap || !primHeap)
42 return 0;
43
44 /* Check if last primary ray has spawned photons (srcIdx >= 0, see
45 * newPhoton()), in which case we write it to the primary heap file
46 * before overwriting it */
47 if (pmap -> lastPrimary.srcIdx >= 0) {
48 if (!fwrite(&pmap -> lastPrimary, sizeof(PhotonPrimary), 1, primHeap))
49 error(SYSTEM, "failed writing photon primary in newPhotonPrimary");
50
51 pmap -> numPrimary++;
52 if (pmap -> numPrimary > PMAP_MAXPRIMARY)
53 error(INTERNAL, "photon primary overflow in newPhotonPrimary");
54 }
55
56 /* Mark unused with negative source index until path spawns a photon (see
57 * newPhoton()) */
58 pmap -> lastPrimary.srcIdx = -1;
59
60 if (primRay) {
61 FVECT dvec;
62
63 /* Reverse incident direction to point to light source */
64 dvec [0] = -primRay -> rdir [0];
65 dvec [1] = -primRay -> rdir [1];
66 dvec [2] = -primRay -> rdir [2];
67 pmap -> lastPrimary.dir = encodedir(dvec);
68 #ifdef PMAP_PRIMARYPOS
69 VCOPY(pmap -> lastPrimary.pos, primRay -> rop);
70 #endif
71 }
72
73 return pmap -> numPrimary;
74 }
75
76
77
78 #ifdef DEBUG_PMAP_CONTRIB
79 static int checkPrimaryHeap (FILE *file)
80 /* Check heap for ordered primaries */
81 {
82 Photon p, lastp;
83 int i, dup;
84
85 rewind(file);
86 memset(&lastp, 0, sizeof(lastp));
87
88 while (fread(&p, sizeof(p), 1, file)) {
89 dup = 1;
90
91 for (i = 0; i <= 2; i++) {
92 if (p.pos [i] < thescene.cuorg [i] ||
93 p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
94
95 sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
96 p.pos [0], p.pos [1], p.pos [2]);
97 error(WARNING, errmsg);
98 }
99
100 dup &= p.pos [i] == lastp.pos [i];
101 }
102
103 if (dup) {
104 sprintf(errmsg,
105 "consecutive duplicate photon in heap at [%f, %f, %f]\n",
106 p.pos [0], p.pos [1], p.pos [2]);
107 error(WARNING, errmsg);
108 }
109 }
110
111 return 0;
112 }
113 #endif
114
115
116
117 static PhotonPrimaryIdx buildPrimaries (PhotonMap *pmap, FILE **primaryHeap,
118 PhotonPrimaryIdx *primaryOfs,
119 unsigned numHeaps)
120 /* Consolidate per-subprocess photon primary heaps into the primary array
121 * pmap -> primaries. Returns offset for primary index linearisation in
122 * numPrimary. The heap files in primaryHeap are closed on return. */
123 {
124 PhotonPrimaryIdx heapLen;
125 unsigned heap;
126
127 if (!pmap || !primaryHeap || !primaryOfs || !numHeaps)
128 return 0;
129
130 pmap -> numPrimary = 0;
131
132 for (heap = 0; heap < numHeaps; heap++) {
133 primaryOfs [heap] = pmap -> numPrimary;
134
135 if (fseek(primaryHeap [heap], 0, SEEK_END))
136 error(SYSTEM, "failed photon primary seek in buildPrimaries");
137 pmap -> numPrimary += heapLen = ftell(primaryHeap [heap]) /
138 sizeof(PhotonPrimary);
139
140 pmap -> primaries = realloc(pmap -> primaries,
141 pmap -> numPrimary *
142 sizeof(PhotonPrimary));
143 if (!pmap -> primaries)
144 error(SYSTEM, "failed photon primary alloc in buildPrimaries");
145
146 rewind(primaryHeap [heap]);
147 if (fread(pmap -> primaries + primaryOfs [heap], sizeof(PhotonPrimary),
148 heapLen, primaryHeap [heap]) != heapLen)
149 error(SYSTEM, "failed reading photon primaries in buildPrimaries");
150
151 fclose(primaryHeap [heap]);
152 }
153
154 return pmap -> numPrimary;
155 }
156
157
158
159 /* Defs for photon emission counter array passed by sub-processes to parent
160 * via shared memory */
161 typedef unsigned long PhotonContribCnt;
162
163 /* Indices for photon emission counter array: num photons stored and num
164 * emitted per source */
165 #define PHOTONCNT_NUMPHOT 0
166 #define PHOTONCNT_NUMEMIT(n) (1 + n)
167
168
169
170 void distribPhotonContrib (PhotonMap* pm, unsigned numProc)
171 {
172 EmissionMap emap;
173 char errmsg2 [128], shmFname [255];
174 unsigned srcIdx, proc;
175 int shmFile, stat, pid;
176 double *srcFlux, /* Emitted flux per light source */
177 srcDistribTarget; /* Target photon count per source */
178 PhotonContribCnt *photonCnt; /* Photon emission counter array */
179 const unsigned photonCntSize = sizeof(PhotonContribCnt) *
180 PHOTONCNT_NUMEMIT(nsources);
181 FILE *primaryHeap [numProc];
182 PhotonPrimaryIdx primaryOfs [numProc];
183
184 if (!pm)
185 error(USER, "no photon map defined in distribPhotonContrib");
186
187 if (!nsources)
188 error(USER, "no light sources in distribPhotonContrib");
189
190 if (nsources > MAXMODLIST)
191 error(USER, "too many light sources in distribPhotonContrib");
192
193 /* Allocate photon flux per light source; this differs for every
194 * source as all sources contribute the same number of distributed
195 * photons (srcDistribTarget), hence the number of photons emitted per
196 * source does not correlate with its emitted flux. The resulting flux
197 * per photon is therefore adjusted individually for each source. */
198 if (!(srcFlux = calloc(nsources, sizeof(double))))
199 error(SYSTEM, "can't allocate source flux in distribPhotonContrib");
200
201 /* ===================================================================
202 * INITIALISATION - Set up emission and scattering funcs
203 * =================================================================== */
204 emap.samples = NULL;
205 emap.src = NULL;
206 emap.maxPartitions = MAXSPART;
207 emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
208 if (!emap.partitions)
209 error(USER, "can't allocate source partitions in distribPhotonContrib");
210
211 /* Initialise contrib photon map */
212 initPhotonMap(pm, PMAP_TYPE_CONTRIB);
213 initPhotonHeap(pm);
214 initPhotonEmissionFuncs();
215 initPhotonScatterFuncs();
216
217 /* Per-subprocess / per-source target counts */
218 pm -> distribTarget /= numProc;
219 srcDistribTarget = nsources ? (double)pm -> distribTarget / nsources : 0;
220
221 /* Get photon ports if specified */
222 if (ambincl == 1)
223 getPhotonPorts();
224
225 /* Get photon sensor modifiers */
226 getPhotonSensors(photonSensorList);
227
228 /* Set up shared mem for photon counters (zeroed by ftruncate) */
229 #if 0
230 snprintf(shmFname, 255, PMAP_SHMFNAME, getpid());
231 shmFile = shm_open(shmFname, O_CREAT | O_RDWR, S_IRUSR | S_IWUSR);
232 #else
233 strcpy(shmFname, PMAP_SHMFNAME);
234 shmFile = mkstemp(shmFname);
235 #endif
236
237 if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0)
238 error(SYSTEM, "failed shared mem init in distribPhotonContrib");
239
240 photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE,
241 MAP_SHARED, shmFile, 0);
242
243 if (photonCnt == MAP_FAILED)
244 error(SYSTEM, "failed shared mem mapping in distribPhotonContrib");
245
246 /* =============================================================
247 * FLUX INTEGRATION - Get total flux emitted from light source
248 * ============================================================= */
249 for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
250 unsigned portCnt = 0;
251
252 srcFlux [srcIdx] = 0;
253 emap.src = source + srcIdx;
254
255 if (photonRepTime)
256 eputs("\n");
257
258 do { /* Need at least one iteration if no ports! */
259 emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
260 : NULL;
261 photonPartition [emap.src -> so -> otype] (&emap);
262
263 if (photonRepTime) {
264 sprintf(errmsg, "Integrating flux from source %s (mod %s) ",
265 source [srcIdx].so -> oname,
266 objptr(source [srcIdx].so -> omod) -> oname);
267
268 if (emap.port) {
269 sprintf(errmsg2, "via port %s ",
270 photonPorts [portCnt].so -> oname);
271 strcat(errmsg, errmsg2);
272 }
273
274 sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
275 strcat(errmsg, errmsg2);
276 eputs(errmsg);
277 fflush(stderr);
278 }
279
280 for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
281 emap.partitionCnt++) {
282 initPhotonEmission(&emap, pdfSamples);
283 srcFlux [srcIdx] += colorAvg(emap.partFlux);
284 }
285
286 portCnt++;
287 } while (portCnt < numPhotonPorts);
288
289 if (srcFlux [srcIdx] < FTINY) {
290 sprintf(errmsg, "source %s has zero emission",
291 source [srcIdx].so -> oname);
292 error(WARNING, errmsg);
293 }
294 }
295
296 if (photonRepTime)
297 eputs("\n");
298
299 /* Init per-subprocess primary heap files */
300 for (proc = 0; proc < numProc; proc++)
301 if (!(primaryHeap [proc] = tmpfile()))
302 error(SYSTEM, "failed opening primary heap file in "
303 "distribPhotonContrib");
304
305 /* MAIN LOOP */
306 for (proc = 0; proc < numProc; proc++) {
307 if (!(pid = fork())) {
308 /* SUBPROCESS ENTERS HERE;
309 * all opened and memory mapped files are inherited */
310
311 /* Local photon counters for this subprocess */
312 unsigned long lastNumPhotons = 0, localNumEmitted = 0;
313 double photonFluxSum = 0; /* Running photon flux sum */
314
315 /* Seed RNGs from PID for decorellated photon distribution */
316 pmapSeed(randSeed + proc, partState);
317 pmapSeed(randSeed + proc, emitState);
318 pmapSeed(randSeed + proc, cntState);
319 pmapSeed(randSeed + proc, mediumState);
320 pmapSeed(randSeed + proc, scatterState);
321 pmapSeed(randSeed + proc, rouletteState);
322
323 /* =============================================================
324 * 2-PASS PHOTON DISTRIBUTION
325 * Pass 1 (pre): emit fraction of target photon count
326 * Pass 2 (main): based on outcome of pass 1, estimate remaining
327 * number of photons to emit to approximate target
328 * count
329 * ============================================================= */
330 for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
331 unsigned portCnt, passCnt = 0, prePassCnt = 0;
332 float srcPreDistrib = preDistrib;
333 double srcNumEmit = 0; /* # to emit from source */
334 unsigned long srcNumDistrib = pm -> numPhotons; /* # stored */
335
336 if (srcFlux [srcIdx] < FTINY)
337 continue;
338
339 while (passCnt < 2) {
340 if (!passCnt) {
341 /* INIT PASS 1 */
342 if (++prePassCnt > maxPreDistrib) {
343 /* Warn if no photons contributed after sufficient
344 * iterations */
345 sprintf(errmsg, "proc %d, source %s: "
346 "too many prepasses, skipped",
347 proc, source [srcIdx].so -> oname);
348 error(WARNING, errmsg);
349 break;
350 }
351
352 /* Num to emit is fraction of target count */
353 srcNumEmit = srcPreDistrib * srcDistribTarget;
354 }
355 else {
356 /* INIT PASS 2 */
357 double srcPhotonFlux, avgPhotonFlux;
358
359 /* Based on the outcome of the predistribution we can now
360 * figure out how many more photons we have to emit from
361 * the current source to meet the target count,
362 * srcDistribTarget. This value is clamped to 0 in case
363 * the target has already been exceeded in pass 1.
364 * srcNumEmit and srcNumDistrib is the number of photons
365 * emitted and distributed (stored) from the current
366 * source in pass 1, respectively. */
367 srcNumDistrib = pm -> numPhotons - srcNumDistrib;
368 srcNumEmit *= srcNumDistrib
369 ? max(srcDistribTarget/srcNumDistrib, 1) - 1
370 : 0;
371
372 if (!srcNumEmit)
373 /* No photons left to distribute in main pass */
374 break;
375
376 srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit;
377 avgPhotonFlux = photonFluxSum / (srcIdx + 1);
378
379 if (avgPhotonFlux > 0 &&
380 srcPhotonFlux / avgPhotonFlux < FTINY) {
381 /* Skip source if its photon flux is grossly below the
382 * running average, indicating negligible contribs at
383 * the expense of excessive distribution time */
384 sprintf(errmsg, "proc %d, source %s: "
385 "itsy bitsy photon flux, skipped",
386 proc, source [srcIdx].so -> oname);
387 error(WARNING, errmsg);
388 srcNumEmit = 0;
389 }
390
391 /* Update sum of photon flux per light source */
392 photonFluxSum += srcPhotonFlux;
393 }
394
395 portCnt = 0;
396 do { /* Need at least one iteration if no ports! */
397 emap.src = source + srcIdx;
398 emap.port = emap.src -> sflags & SDISTANT
399 ? photonPorts + portCnt : NULL;
400 photonPartition [emap.src -> so -> otype] (&emap);
401
402 if (photonRepTime && !proc) {
403 if (!passCnt)
404 sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
405 prePassCnt, source [srcIdx].so -> oname,
406 objptr(source[srcIdx].so->omod) -> oname);
407 else
408 sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
409 source [srcIdx].so -> oname,
410 objptr(source[srcIdx].so->omod) -> oname);
411
412 if (emap.port) {
413 sprintf(errmsg2, "via port %s ",
414 photonPorts [portCnt].so -> oname);
415 strcat(errmsg, errmsg2);
416 }
417
418 sprintf(errmsg2, "(%lu partitions)\n",
419 emap.numPartitions);
420 strcat(errmsg, errmsg2);
421 eputs(errmsg);
422 fflush(stderr);
423 }
424
425 for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
426 emap.partitionCnt++) {
427 double partNumEmit;
428 unsigned long partEmitCnt;
429
430 /* Get photon origin within current source partishunn
431 * and build emission map */
432 photonOrigin [emap.src -> so -> otype] (&emap);
433 initPhotonEmission(&emap, pdfSamples);
434
435 /* Number of photons to emit from ziss partishunn;
436 * scale according to its normalised contribushunn to
437 * the emitted source flux */
438 partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
439 srcFlux [srcIdx];
440 partEmitCnt = (unsigned long)partNumEmit;
441
442 /* Probabilistically account for fractional photons */
443 if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
444 partEmitCnt++;
445
446 /* Update local and shared global emission counter */
447 localNumEmitted += partEmitCnt;
448 photonCnt [PHOTONCNT_NUMEMIT(srcIdx)] += partEmitCnt;
449
450 /* Integer counter avoids FP rounding errors */
451 while (partEmitCnt--) {
452 RAY photonRay;
453
454 /* Emit photon according to PDF (if any), allocate
455 * associated primary ray, and trace through scene
456 * until absorbed/leaked; emitPhoton() sets the
457 * emitting light source index in photonRay */
458 emitPhoton(&emap, &photonRay);
459 newPhotonPrimary(pm, &photonRay, primaryHeap[proc]);
460 /* Set subprocess index in photonRay for post-
461 * distrib primary index linearisation; this is
462 * propagated with the primary index in photonRay
463 * and set for photon hits by newPhoton() */
464 PMAP_SETRAYPROC(&photonRay, proc);
465 tracePhoton(&photonRay);
466 }
467
468 /* Update shared global photon count */
469 photonCnt [PHOTONCNT_NUMPHOT] += pm -> numPhotons -
470 lastNumPhotons;
471 lastNumPhotons = pm -> numPhotons;
472 }
473
474 portCnt++;
475 } while (portCnt < numPhotonPorts);
476
477 if (pm -> numPhotons == srcNumDistrib)
478 /* Double predistrib factor in case no photons were stored
479 * for this source and redo pass 1 */
480 srcPreDistrib *= 2;
481 else {
482 /* Now do pass 2 */
483 passCnt++;
484 /* if (photonRepTime)
485 eputs("\n"); */
486 }
487 }
488 }
489
490 /* Flush heap buffa one final time to prevent data corruption */
491 flushPhotonHeap(pm);
492 fclose(pm -> heap);
493
494 /* Flush final photon primary to primary heap file */
495 newPhotonPrimary(pm, NULL, primaryHeap [proc]);
496 fclose(primaryHeap [proc]);
497
498 #ifdef DEBUG_PMAP
499 sprintf(errmsg, "Proc %d exited with total %ld photons\n", proc,
500 pm -> numPhotons);
501 eputs(errmsg);
502 #endif
503
504 exit(0);
505 }
506 else if (pid < 0)
507 error(SYSTEM, "failed to fork subprocess in distribPhotonContrib");
508 }
509
510 /* PARENT PROCESS CONTINUES HERE */
511 /* Record start time and enable progress report signal handler */
512 repStartTime = time(NULL);
513 #ifdef SIGCONT
514 signal(SIGCONT, pmapDistribReport);
515 #endif
516 /*
517 if (photonRepTime)
518 eputs("\n"); */
519
520 /* Wait for subprocesses to complete while reporting progress */
521 proc = numProc;
522 while (proc) {
523 while (waitpid(-1, &stat, WNOHANG) > 0) {
524 /* Subprocess exited; check status */
525 if (!WIFEXITED(stat) || WEXITSTATUS(stat))
526 error(USER, "failed photon distribution");
527
528 --proc;
529 }
530
531 /* Nod off for a bit and update progress */
532 sleep(1);
533
534 /* Update progress report from shared subprocess counters */
535 repComplete = pm -> distribTarget * numProc;
536 repProgress = photonCnt [PHOTONCNT_NUMPHOT];
537 for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++)
538 repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
539
540 /* Get global photon count from shmem updated by subprocs */
541 pm -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT];
542
543 if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
544 pmapDistribReport();
545 #ifdef SIGCONT
546 else signal(SIGCONT, pmapDistribReport);
547 #endif
548 }
549
550 /* ================================================================
551 * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
552 * ================================================================ */
553 #ifdef SIGCONT
554 signal(SIGCONT, SIG_DFL);
555 #endif
556 free(emap.samples);
557
558 if (!pm -> numPhotons)
559 error(USER, "empty photon map");
560
561 /* Load per-subprocess primary rays into pm -> primary array */
562 pm -> numPrimary = buildPrimaries(pm, primaryHeap, primaryOfs, numProc);
563 if (!pm -> numPrimary)
564 error(INTERNAL, "no primary rays in contribution photon map");
565
566 /* Set photon flux per source */
567 for (srcIdx = 0; srcIdx < nsources; srcIdx++)
568 srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
569
570 /* Photon counters no longer needed, unmap shared memory */
571 munmap(photonCnt, sizeof(*photonCnt));
572 close(shmFile);
573 #if 0
574 shm_unlink(shmFname);
575 #else
576 unlink(shmFname);
577 #endif
578
579 if (photonRepTime) {
580 eputs("\nBuilding contrib photon map...\n");
581 fflush(stderr);
582 }
583
584 /* Build underlying data structure; heap is destroyed */
585 buildPhotonMap(pm, srcFlux, primaryOfs, numProc);
586 }