ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmap.c
Revision: 2.13
Committed: Mon Aug 14 21:12:10 2017 UTC (6 years, 9 months ago) by rschregle
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R1
Changes since 2.12: +170 -102 lines
Log Message:
Updated photon map code for Windows; no multproc or ooC for now

File Contents

# User Rev Content
1 greg 2.9 #ifndef lint
2 rschregle 2.13 static const char RCSid[] = "$Id: pmap.c,v 2.12 2016/09/26 20:19:30 greg Exp $";
3 greg 2.9 #endif
4 rschregle 2.11
5 rschregle 2.13
6 greg 2.1 /*
7 rschregle 2.11 ======================================================================
8 greg 2.1 Photon map main module
9    
10     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
11     (c) Fraunhofer Institute for Solar Energy Systems,
12 rschregle 2.4 (c) Lucerne University of Applied Sciences and Arts,
13 rschregle 2.11 supported by the Swiss National Science Foundation (SNSF, #147053)
14     ======================================================================
15 greg 2.1
16 rschregle 2.13 $Id: pmap.c,v 2.12 2016/09/26 20:19:30 greg Exp $
17 greg 2.1 */
18    
19    
20     #include "pmap.h"
21     #include "pmapmat.h"
22     #include "pmapsrc.h"
23     #include "pmaprand.h"
24     #include "pmapio.h"
25     #include "pmapbias.h"
26     #include "pmapdiag.h"
27     #include "otypes.h"
28     #include <time.h>
29 rschregle 2.13 #if NIX
30     #include <sys/stat.h>
31     #include <sys/mman.h>
32     #include <sys/wait.h>
33     #endif
34 greg 2.1
35    
36     void savePmaps (const PhotonMap **pmaps, int argc, char **argv)
37     {
38     unsigned t;
39    
40     for (t = 0; t < NUM_PMAP_TYPES; t++) {
41     if (pmaps [t])
42 greg 2.7 savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv);
43 greg 2.1 }
44     }
45    
46    
47    
48     static int photonParticipate (RAY *ray)
49     /* Trace photon through participating medium. Returns 1 if passed through,
50     or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */
51     {
52     int i;
53     RREAL cosTheta, cosPhi, du, dv;
54     const float cext = colorAvg(ray -> cext),
55     albedo = colorAvg(ray -> albedo);
56     FVECT u, v;
57     COLOR cvext;
58    
59     /* Mean free distance until interaction with medium */
60     ray -> rmax = -log(pmapRandom(mediumState)) / cext;
61    
62     while (!localhit(ray, &thescene)) {
63     setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]),
64     exp(-ray -> rmax * ray -> cext [1]),
65     exp(-ray -> rmax * ray -> cext [2]));
66    
67     /* Modify ray color and normalise */
68     multcolor(ray -> rcol, cvext);
69     colorNorm(ray -> rcol);
70     VCOPY(ray -> rorg, ray -> rop);
71    
72 rschregle 2.11 if (albedo > FTINY && ray -> rlvl > 0)
73 greg 2.1 /* Add to volume photon map */
74 rschregle 2.11 newPhoton(volumePmap, ray);
75 greg 2.1
76     /* Absorbed? */
77 rschregle 2.11 if (pmapRandom(rouletteState) > albedo)
78     return 0;
79 greg 2.1
80     /* Colour bleeding without attenuation (?) */
81     multcolor(ray -> rcol, ray -> albedo);
82     scalecolor(ray -> rcol, 1 / albedo);
83    
84     /* Scatter photon */
85     cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1
86     : 1 / (2 * ray -> gecc) *
87     (1 + ray -> gecc * ray -> gecc -
88     (1 - ray -> gecc * ray -> gecc) /
89     (1 - ray -> gecc + 2 * ray -> gecc *
90     pmapRandom(scatterState)));
91    
92     cosPhi = cos(2 * PI * pmapRandom(scatterState));
93     du = dv = sqrt(1 - cosTheta * cosTheta); /* sin(theta) */
94     du *= cosPhi;
95     dv *= sqrt(1 - cosPhi * cosPhi); /* sin(phi) */
96    
97     /* Get axes u & v perpendicular to photon direction */
98     i = 0;
99     do {
100     v [0] = v [1] = v [2] = 0;
101     v [i++] = 1;
102     fcross(u, v, ray -> rdir);
103     } while (normalize(u) < FTINY);
104     fcross(v, ray -> rdir, u);
105    
106     for (i = 0; i < 3; i++)
107     ray -> rdir [i] = du * u [i] + dv * v [i] +
108     cosTheta * ray -> rdir [i];
109     ray -> rlvl++;
110     ray -> rmax = -log(pmapRandom(mediumState)) / cext;
111     }
112    
113     setcolor(cvext, exp(-ray -> rot * ray -> cext [0]),
114     exp(-ray -> rot * ray -> cext [1]),
115     exp(-ray -> rot * ray -> cext [2]));
116    
117     /* Modify ray color and normalise */
118     multcolor(ray -> rcol, cvext);
119     colorNorm(ray -> rcol);
120    
121     /* Passed through medium */
122     return 1;
123     }
124    
125    
126    
127     void tracePhoton (RAY *ray)
128     /* Follow photon as it bounces around... */
129     {
130     long mod;
131 rschregle 2.13 OBJREC *mat, *port = NULL;
132    
133     if (!ray -> parent) {
134     /* !!! PHOTON PORT REJECTION SAMPLING HACK: get photon port for
135     * !!! primary ray from ray -> ro, then reset the latter to NULL so
136     * !!! as not to interfere with localhit() */
137     port = ray -> ro;
138     ray -> ro = NULL;
139     }
140 greg 2.1
141     if (ray -> rlvl > photonMaxBounce) {
142 rschregle 2.5 #ifdef PMAP_RUNAWAY_WARN
143 greg 2.1 error(WARNING, "runaway photon!");
144 rschregle 2.5 #endif
145 greg 2.1 return;
146     }
147 rschregle 2.5
148 greg 2.1 if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
149     return;
150 rschregle 2.13
151 greg 2.1 if (localhit(ray, &thescene)) {
152     mod = ray -> ro -> omod;
153 rschregle 2.13
154     if (port && ray -> ro != port) {
155     /* !!! PHOTON PORT REJECTION SAMPLING HACK !!!
156     * Terminate photon if emitted from port without intersecting it;
157     * this can happen when the port's partitions extend beyond its
158     * actual geometry, e.g. with polygons. Since the total flux
159     * relayed by the port is based on the (in this case) larger
160     * partition area, it is overestimated; terminating these photons
161     * constitutes rejection sampling and thereby compensates any bias
162     * incurred by the overestimated flux. */
163     #ifdef PMAP_PORTREJECT_WARN
164     sprintf(errmsg, "photon outside port %s", ray -> ro -> oname);
165     error(WARNING, errmsg);
166     #endif
167     return;
168     }
169    
170 greg 2.1 if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) {
171     /* Transfer ray if modifier is VOID or clipped within antimatta */
172     RAY tray;
173     photonRay(ray, &tray, PMAP_XFER, NULL);
174     tracePhoton(&tray);
175     }
176     else {
177     /* Scatter for modifier material */
178     mat = objptr(mod);
179     photonScatter [mat -> otype] (mat, ray);
180     }
181     }
182     }
183    
184    
185    
186     static void preComputeGlobal (PhotonMap *pmap)
187 rschregle 2.11 /* Precompute irradiance from global photons for final gathering for
188     a random subset of finalGather * pmap -> numPhotons photons, and builds
189     the photon map, discarding the original photons. */
190     /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */
191     {
192     unsigned long i, numPreComp;
193     unsigned j;
194     PhotonIdx pIdx;
195     Photon photon;
196     RAY ray;
197     PhotonMap nuPmap;
198 greg 2.1
199 rschregle 2.11 repComplete = numPreComp = finalGather * pmap -> numPhotons;
200 greg 2.1
201 rschregle 2.13 if (verbose) {
202     sprintf(errmsg,
203     "\nPrecomputing irradiance for %ld global photons\n",
204 rschregle 2.11 numPreComp);
205 greg 2.1 eputs(errmsg);
206 rschregle 2.13 #if NIX
207 greg 2.1 fflush(stderr);
208 rschregle 2.13 #endif
209 greg 2.1 }
210    
211 rschregle 2.11 /* Copy photon map for precomputed photons */
212     memcpy(&nuPmap, pmap, sizeof(PhotonMap));
213    
214     /* Zero counters, init new heap and extents */
215     nuPmap.numPhotons = 0;
216     initPhotonHeap(&nuPmap);
217    
218     for (j = 0; j < 3; j++) {
219     nuPmap.minPos [j] = FHUGE;
220     nuPmap.maxPos [j] = -FHUGE;
221 greg 2.1 }
222 rschregle 2.11
223 greg 2.1 /* Record start time, baby */
224     repStartTime = time(NULL);
225 rschregle 2.11 #ifdef SIGCONT
226     signal(SIGCONT, pmapPreCompReport);
227     #endif
228 greg 2.1 repProgress = 0;
229    
230 rschregle 2.11 photonRay(NULL, &ray, PRIMARY, NULL);
231     ray.ro = NULL;
232    
233     for (i = 0; i < numPreComp; i++) {
234     /* Get random photon from stratified distribution in source heap to
235 rschregle 2.13 * avoid duplicates and clustering */
236 rschregle 2.11 pIdx = firstPhoton(pmap) +
237     (unsigned long)((i + pmapRandom(pmap -> randState)) /
238     finalGather);
239     getPhoton(pmap, pIdx, &photon);
240    
241     /* Init dummy photon ray with intersection at photon position */
242     VCOPY(ray.rop, photon.pos);
243     for (j = 0; j < 3; j++)
244     ray.ron [j] = photon.norm [j] / 127.0;
245    
246     /* Get density estimate at photon position */
247     photonDensity(pmap, &ray, ray.rcol);
248    
249     /* Append photon to new heap from ray */
250     newPhoton(&nuPmap, &ray);
251 greg 2.1
252 rschregle 2.11 /* Update progress */
253 greg 2.1 repProgress++;
254    
255     if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
256     pmapPreCompReport();
257 rschregle 2.11 #ifdef SIGCONT
258     else signal(SIGCONT, pmapPreCompReport);
259     #endif
260 greg 2.1 }
261    
262 rschregle 2.11 /* Flush heap */
263     flushPhotonHeap(&nuPmap);
264    
265     #ifdef SIGCONT
266     signal(SIGCONT, SIG_DFL);
267     #endif
268    
269     /* Trash original pmap, replace with precomputed one */
270     deletePhotons(pmap);
271     memcpy(pmap, &nuPmap, sizeof(PhotonMap));
272 greg 2.1
273 rschregle 2.13 if (verbose) {
274     eputs("\nRebuilding precomputed photon map\n");
275     #if NIX
276 greg 2.1 fflush(stderr);
277 rschregle 2.13 #endif
278 greg 2.1 }
279 rschregle 2.11
280     /* Rebuild underlying data structure, destroying heap */
281     buildPhotonMap(pmap, NULL, NULL, 1);
282 greg 2.1 }
283    
284    
285    
286 rschregle 2.11 typedef struct {
287     unsigned long numPhotons [NUM_PMAP_TYPES],
288     numEmitted, numComplete;
289     } PhotonCnt;
290    
291    
292    
293     void distribPhotons (PhotonMap **pmaps, unsigned numProc)
294     {
295     EmissionMap emap;
296 rschregle 2.13 char errmsg2 [128], shmFname [PMAP_TMPFNLEN];
297 rschregle 2.11 unsigned t, srcIdx, proc;
298     double totalFlux = 0;
299     int shmFile, stat, pid;
300     PhotonMap *pm;
301     PhotonCnt *photonCnt;
302 greg 2.1
303 rschregle 2.8 for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++);
304 rschregle 2.11
305 greg 2.1 if (t >= NUM_PMAP_TYPES)
306 rschregle 2.11 error(USER, "no photon maps defined in distribPhotons");
307 greg 2.1
308     if (!nsources)
309 rschregle 2.11 error(USER, "no light sources in distribPhotons");
310 greg 2.1
311     /* ===================================================================
312     * INITIALISATION - Set up emission and scattering funcs
313     * =================================================================== */
314     emap.samples = NULL;
315     emap.maxPartitions = MAXSPART;
316     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
317     if (!emap.partitions)
318 rschregle 2.11 error(INTERNAL, "can't allocate source partitions in distribPhotons");
319 greg 2.1
320     /* Initialise all defined photon maps */
321     for (t = 0; t < NUM_PMAP_TYPES; t++)
322 rschregle 2.11 if (pmaps [t]) {
323     initPhotonMap(pmaps [t], t);
324     /* Open photon heapfile */
325     initPhotonHeap(pmaps [t]);
326     /* Per-subprocess target count */
327     pmaps [t] -> distribTarget /= numProc;
328 rschregle 2.13
329     if (!pmaps [t] -> distribTarget)
330     error(INTERNAL, "no photons to distribute in distribPhotons");
331 rschregle 2.11 }
332 greg 2.1
333     initPhotonEmissionFuncs();
334     initPhotonScatterFuncs();
335    
336     /* Get photon ports if specified */
337     if (ambincl == 1)
338     getPhotonPorts();
339    
340     /* Get photon sensor modifiers */
341     getPhotonSensors(photonSensorList);
342    
343 rschregle 2.13 #if NIX
344 rschregle 2.11 /* Set up shared mem for photon counters (zeroed by ftruncate) */
345 rschregle 2.13 strcpy(shmFname, PMAP_TMPFNAME);
346 rschregle 2.11 shmFile = mkstemp(shmFname);
347    
348 rschregle 2.13 if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0)
349     error(SYSTEM, "failed shared mem init in distribPhotons");
350 rschregle 2.11
351     photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE,
352     MAP_SHARED, shmFile, 0);
353    
354     if (photonCnt == MAP_FAILED)
355 rschregle 2.13 error(SYSTEM, "failed mapping shared memory in distribPhotons");
356     #else
357     /* Allocate photon counters statically on Windoze */
358     if (!(photonCnt = malloc(sizeof(PhotonCnt))))
359     error(SYSTEM, "failed trivial malloc in distribPhotons");
360     photonCnt -> numEmitted = photonCnt -> numComplete = 0;
361     #endif /* NIX */
362    
363     if (verbose) {
364     sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
365    
366     if (photonPorts) {
367     sprintf(errmsg2, " via %d ports", numPhotonPorts);
368     strcat(errmsg, errmsg2);
369     }
370    
371     strcat(errmsg, "\n");
372     eputs(errmsg);
373     }
374 greg 2.1
375     /* ===================================================================
376     * FLUX INTEGRATION - Get total photon flux from light sources
377     * =================================================================== */
378 rschregle 2.11 for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
379 greg 2.1 unsigned portCnt = 0;
380     emap.src = source + srcIdx;
381    
382 rschregle 2.11 do { /* Need at least one iteration if no ports! */
383 greg 2.1 emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
384     : NULL;
385     photonPartition [emap.src -> so -> otype] (&emap);
386    
387 rschregle 2.13 if (verbose) {
388     sprintf(errmsg, "\tIntegrating flux from source %s ",
389 greg 2.1 source [srcIdx].so -> oname);
390 rschregle 2.13
391 greg 2.1 if (emap.port) {
392     sprintf(errmsg2, "via port %s ",
393     photonPorts [portCnt].so -> oname);
394     strcat(errmsg, errmsg2);
395     }
396 rschregle 2.13
397     sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
398 greg 2.1 strcat(errmsg, errmsg2);
399     eputs(errmsg);
400 rschregle 2.13 #if NIX
401 greg 2.1 fflush(stderr);
402 rschregle 2.13 #endif
403 greg 2.1 }
404    
405     for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
406     emap.partitionCnt++) {
407     initPhotonEmission(&emap, pdfSamples);
408     totalFlux += colorAvg(emap.partFlux);
409     }
410    
411     portCnt++;
412     } while (portCnt < numPhotonPorts);
413     }
414    
415     if (totalFlux < FTINY)
416     error(USER, "zero flux from light sources");
417 rschregle 2.13
418     /* Record start time for progress reports */
419     repStartTime = time(NULL);
420    
421     if (verbose) {
422     sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
423     eputs(errmsg);
424     }
425 greg 2.1
426 rschregle 2.11 /* MAIN LOOP */
427     for (proc = 0; proc < numProc; proc++) {
428 rschregle 2.13 #if NIX
429 rschregle 2.11 if (!(pid = fork())) {
430 rschregle 2.13 /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */
431     #else
432     if (1) {
433     /* No subprocess under Windoze */
434     #endif
435     /* Local photon counters for this subprocess */
436 rschregle 2.11 unsigned passCnt = 0, prePassCnt = 0;
437     unsigned long lastNumPhotons [NUM_PMAP_TYPES];
438     unsigned long localNumEmitted = 0; /* Num photons emitted by this
439     subprocess alone */
440 greg 2.1
441 rschregle 2.11 /* Seed RNGs from PID for decorellated photon distribution */
442     pmapSeed(randSeed + proc, partState);
443 rschregle 2.13 pmapSeed(randSeed + (proc + 1) % numProc, emitState);
444     pmapSeed(randSeed + (proc + 2) % numProc, cntState);
445     pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
446     pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
447     pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
448 greg 2.1
449     for (t = 0; t < NUM_PMAP_TYPES; t++)
450 rschregle 2.11 lastNumPhotons [t] = 0;
451    
452     /* =============================================================
453     * 2-PASS PHOTON DISTRIBUTION
454     * Pass 1 (pre): emit fraction of target photon count
455     * Pass 2 (main): based on outcome of pass 1, estimate remaining
456     * number of photons to emit to approximate target
457     * count
458     * ============================================================= */
459 greg 2.1 do {
460 rschregle 2.11 double numEmit;
461 greg 2.1
462 rschregle 2.11 if (!passCnt) {
463     /* INIT PASS 1 */
464     /* Skip if no photons contributed after sufficient
465     * iterations; make it clear to user which photon maps are
466     * missing so (s)he can check geometry and materials */
467     if (++prePassCnt > maxPreDistrib) {
468 rschregle 2.13 sprintf(errmsg, "proc %d: too many prepasses", proc);
469 rschregle 2.11
470     for (t = 0; t < NUM_PMAP_TYPES; t++)
471     if (pmaps [t] && !pmaps [t] -> numPhotons) {
472     sprintf(errmsg2, ", no %s photons stored",
473     pmapName [t]);
474     strcat(errmsg, errmsg2);
475     }
476    
477     error(USER, errmsg);
478     break;
479 greg 2.1 }
480 rschregle 2.11
481     /* Num to emit is fraction of minimum target count */
482     numEmit = FHUGE;
483 greg 2.1
484 rschregle 2.11 for (t = 0; t < NUM_PMAP_TYPES; t++)
485     if (pmaps [t])
486     numEmit = min(pmaps [t] -> distribTarget, numEmit);
487    
488     numEmit *= preDistrib;
489 greg 2.1 }
490 rschregle 2.11 else {
491     /* INIT PASS 2 */
492     /* Based on the outcome of the predistribution we can now
493     * estimate how many more photons we have to emit for each
494     * photon map to meet its respective target count. This
495     * value is clamped to 0 in case the target has already been
496     * exceeded in the pass 1. */
497     double maxDistribRatio = 0;
498    
499     /* Set the distribution ratio for each map; this indicates
500     * how many photons of each respective type are stored per
501     * emitted photon, and is used as probability for storing a
502     * photon by newPhoton(). Since this biases the photon
503     * density, newPhoton() promotes the flux of stored photons
504     * to compensate. */
505     for (t = 0; t < NUM_PMAP_TYPES; t++)
506     if ((pm = pmaps [t])) {
507     pm -> distribRatio = (double)pm -> distribTarget /
508     pm -> numPhotons - 1;
509    
510     /* Check if photon map "overflowed", i.e. exceeded its
511     * target count in the prepass; correcting the photon
512     * flux via the distribution ratio is no longer
513     * possible, as no more photons of this type will be
514     * stored, so notify the user rather than deliver
515     * incorrect results. In future we should handle this
516     * more intelligently by using the photonFlux in each
517     * photon map to individually correct the flux after
518     * distribution. */
519     if (pm -> distribRatio <= FTINY) {
520     sprintf(errmsg, "%s photon map overflow in "
521     "prepass, reduce -apD", pmapName [t]);
522     error(INTERNAL, errmsg);
523     }
524    
525     maxDistribRatio = max(pm -> distribRatio,
526     maxDistribRatio);
527     }
528 greg 2.1
529 rschregle 2.11 /* Normalise distribution ratios and calculate number of
530     * photons to emit in main pass */
531     for (t = 0; t < NUM_PMAP_TYPES; t++)
532     if ((pm = pmaps [t]))
533     pm -> distribRatio /= maxDistribRatio;
534    
535     if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY)
536     /* No photons left to distribute in main pass */
537     break;
538     }
539    
540 rschregle 2.13 /* Update shared completion counter for progreport by parent */
541 rschregle 2.11 photonCnt -> numComplete += numEmit;
542    
543     /* PHOTON DISTRIBUTION LOOP */
544     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
545     unsigned portCnt = 0;
546     emap.src = source + srcIdx;
547    
548     do { /* Need at least one iteration if no ports! */
549     emap.port = emap.src -> sflags & SDISTANT
550     ? photonPorts + portCnt : NULL;
551     photonPartition [emap.src -> so -> otype] (&emap);
552    
553 rschregle 2.13 if (verbose && !proc) {
554     /* Output from subproc 0 only to avoid race condition
555     * on console I/O */
556 rschregle 2.11 if (!passCnt)
557 rschregle 2.13 sprintf(errmsg, "\tPREPASS %d on source %s ",
558 rschregle 2.11 prePassCnt, source [srcIdx].so -> oname);
559     else
560 rschregle 2.13 sprintf(errmsg, "\tMAIN PASS on source %s ",
561 rschregle 2.11 source [srcIdx].so -> oname);
562 rschregle 2.13
563 rschregle 2.11 if (emap.port) {
564     sprintf(errmsg2, "via port %s ",
565     photonPorts [portCnt].so -> oname);
566     strcat(errmsg, errmsg2);
567     }
568 rschregle 2.13
569 rschregle 2.11 sprintf(errmsg2, "(%lu partitions)\n",
570     emap.numPartitions);
571     strcat(errmsg, errmsg2);
572     eputs(errmsg);
573 rschregle 2.13 #if NIX
574 rschregle 2.11 fflush(stderr);
575 rschregle 2.13 #endif
576 rschregle 2.11 }
577 greg 2.1
578 rschregle 2.11 for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
579     emap.partitionCnt++) {
580     double partNumEmit;
581     unsigned long partEmitCnt;
582    
583     /* Get photon origin within current source partishunn
584     * and build emission map */
585     photonOrigin [emap.src -> so -> otype] (&emap);
586     initPhotonEmission(&emap, pdfSamples);
587    
588     /* Number of photons to emit from ziss partishunn --
589     * proportional to flux; photon ray weight and scalar
590 rschregle 2.13 * flux are uniform (latter only varying in RGB). */
591 rschregle 2.11 partNumEmit = numEmit * colorAvg(emap.partFlux) /
592     totalFlux;
593     partEmitCnt = (unsigned long)partNumEmit;
594    
595     /* Probabilistically account for fractional photons */
596     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
597     partEmitCnt++;
598    
599     /* Update local and shared (global) emission counter */
600     photonCnt -> numEmitted += partEmitCnt;
601     localNumEmitted += partEmitCnt;
602    
603     /* Integer counter avoids FP rounding errors during
604     * iteration */
605     while (partEmitCnt--) {
606     RAY photonRay;
607    
608     /* Emit photon based on PDF and trace through scene
609     * until absorbed/leaked */
610     emitPhoton(&emap, &photonRay);
611 rschregle 2.13 #if 1
612     if (emap.port)
613     /* !!! PHOTON PORT REJECTION SAMPLING HACK: set
614     * !!! photon port as fake hit object for
615     * !!! primary ray to check for intersection in
616     * !!! tracePhoton() */
617     photonRay.ro = emap.port -> so;
618     #endif
619 rschregle 2.11 tracePhoton(&photonRay);
620     }
621 rschregle 2.13
622 rschregle 2.11 /* Update shared global photon count for each pmap */
623     for (t = 0; t < NUM_PMAP_TYPES; t++)
624     if (pmaps [t]) {
625     photonCnt -> numPhotons [t] +=
626     pmaps [t] -> numPhotons - lastNumPhotons [t];
627     lastNumPhotons [t] = pmaps [t] -> numPhotons;
628     }
629 rschregle 2.13 #if !NIX
630     /* Synchronous progress report on Windoze */
631     if (!proc && photonRepTime > 0 &&
632     time(NULL) >= repLastTime + photonRepTime) {
633     repEmitted = repProgress = photonCnt -> numEmitted;
634     repComplete = photonCnt -> numComplete;
635     pmapDistribReport();
636     }
637     #endif
638 rschregle 2.11 }
639 greg 2.1
640 rschregle 2.11 portCnt++;
641     } while (portCnt < numPhotonPorts);
642     }
643    
644     for (t = 0; t < NUM_PMAP_TYPES; t++)
645     if (pmaps [t] && !pmaps [t] -> numPhotons) {
646     /* Double preDistrib in case a photon map is empty and
647     * redo pass 1 --> possibility of infinite loop for
648     * pathological scenes (e.g. absorbing materials) */
649     preDistrib *= 2;
650     break;
651 greg 2.1 }
652 rschregle 2.11
653 rschregle 2.13 if (t >= NUM_PMAP_TYPES)
654 rschregle 2.11 /* No empty photon maps found; now do pass 2 */
655     passCnt++;
656     } while (passCnt < 2);
657    
658 rschregle 2.13 /* Flush heap buffa for every pmap one final time;
659     * avoids potential data corruption! */
660 rschregle 2.11 for (t = 0; t < NUM_PMAP_TYPES; t++)
661     if (pmaps [t]) {
662     flushPhotonHeap(pmaps [t]);
663 rschregle 2.13 /* Heap file closed automatically on exit
664     fclose(pmaps [t] -> heap); */
665 rschregle 2.11 #ifdef DEBUG_PMAP
666 rschregle 2.13 sprintf(errmsg, "Proc %d: total %ld photons\n", proc,
667 rschregle 2.11 pmaps [t] -> numPhotons);
668     eputs(errmsg);
669     #endif
670     }
671 rschregle 2.13 #if NIX
672     /* Terminate subprocess */
673 rschregle 2.11 exit(0);
674 rschregle 2.13 #endif
675 greg 2.1 }
676 rschregle 2.11 else if (pid < 0)
677     error(SYSTEM, "failed to fork subprocess in distribPhotons");
678     }
679    
680 rschregle 2.13 #if NIX
681 rschregle 2.11 /* PARENT PROCESS CONTINUES HERE */
682     #ifdef SIGCONT
683 rschregle 2.13 /* Enable progress report signal handler */
684 rschregle 2.11 signal(SIGCONT, pmapDistribReport);
685 rschregle 2.13 #endif
686     /* Wait for subprocesses complete while reporting progress */
687 rschregle 2.11 proc = numProc;
688     while (proc) {
689     while (waitpid(-1, &stat, WNOHANG) > 0) {
690     /* Subprocess exited; check status */
691     if (!WIFEXITED(stat) || WEXITSTATUS(stat))
692     error(USER, "failed photon distribution");
693    
694     --proc;
695     }
696    
697     /* Nod off for a bit and update progress */
698     sleep(1);
699 rschregle 2.13
700     /* Asynchronous progress report from shared subprocess counters */
701 rschregle 2.11 repEmitted = repProgress = photonCnt -> numEmitted;
702 rschregle 2.13 repComplete = photonCnt -> numComplete;
703 rschregle 2.11
704 rschregle 2.13 repProgress = repComplete = 0;
705 greg 2.1 for (t = 0; t < NUM_PMAP_TYPES; t++)
706 rschregle 2.11 if ((pm = pmaps [t])) {
707     /* Get global photon count from shmem updated by subprocs */
708 rschregle 2.13 repProgress += pm -> numPhotons = photonCnt -> numPhotons [t];
709     repComplete += pm -> distribTarget;
710 greg 2.1 }
711 rschregle 2.13 repComplete *= numProc;
712 rschregle 2.11
713     if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
714     pmapDistribReport();
715     #ifdef SIGCONT
716     else signal(SIGCONT, pmapDistribReport);
717     #endif
718     }
719 rschregle 2.13 #endif /* NIX */
720 greg 2.1
721     /* ===================================================================
722 rschregle 2.11 * POST-DISTRIBUTION - Set photon flux and build data struct for photon
723     * storage, etc.
724 greg 2.1 * =================================================================== */
725 rschregle 2.11 #ifdef SIGCONT
726 rschregle 2.13 /* Reset signal handler */
727 rschregle 2.11 signal(SIGCONT, SIG_DFL);
728     #endif
729 greg 2.1 free(emap.samples);
730    
731 rschregle 2.13 /* Set photon flux */
732 rschregle 2.11 totalFlux /= photonCnt -> numEmitted;
733 rschregle 2.13 #if NIX
734 rschregle 2.11 /* Photon counters no longer needed, unmap shared memory */
735     munmap(photonCnt, sizeof(*photonCnt));
736     close(shmFile);
737 rschregle 2.13 unlink(shmFname);
738 rschregle 2.11 #else
739 rschregle 2.13 free(photonCnt);
740 rschregle 2.11 #endif
741 rschregle 2.13 if (verbose)
742     eputs("\n");
743    
744 greg 2.1 for (t = 0; t < NUM_PMAP_TYPES; t++)
745 rschregle 2.8 if (pmaps [t]) {
746 rschregle 2.13 if (verbose) {
747     sprintf(errmsg, "Building %s photon map\n", pmapName [t]);
748 greg 2.1 eputs(errmsg);
749 rschregle 2.13 #if NIX
750 greg 2.1 fflush(stderr);
751 rschregle 2.13 #endif
752 greg 2.1 }
753 rschregle 2.11
754     /* Build underlying data structure; heap is destroyed */
755     buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
756 greg 2.1 }
757 rschregle 2.13
758 greg 2.1 /* Precompute photon irradiance if necessary */
759 rschregle 2.13 if (preCompPmap) {
760     if (verbose)
761     eputs("\n");
762 greg 2.1 preComputeGlobal(preCompPmap);
763 rschregle 2.13 }
764    
765     if (verbose)
766     eputs("\n");
767 greg 2.1 }