ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmap.c
Revision: 2.14
Committed: Tue Mar 20 19:55:33 2018 UTC (6 years, 1 month ago) by rschregle
Content type: text/plain
Branch: MAIN
Changes since 2.13: +4 -5 lines
Log Message:
Added -ae/-ai ambient exclude options to mkpmap, cleaned up opt parsing.

File Contents

# User Rev Content
1 greg 2.9 #ifndef lint
2 rschregle 2.14 static const char RCSid[] = "$Id: pmap.c,v 2.13 2017/08/14 21:12:10 rschregle 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.14 $Id: pmap.c,v 2.13 2017/08/14 21:12:10 rschregle 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 rschregle 2.14 /* Get photon ports from modifier list */
337     getPhotonPorts(photonPortList);
338 greg 2.1
339     /* Get photon sensor modifiers */
340     getPhotonSensors(photonSensorList);
341    
342 rschregle 2.13 #if NIX
343 rschregle 2.11 /* Set up shared mem for photon counters (zeroed by ftruncate) */
344 rschregle 2.13 strcpy(shmFname, PMAP_TMPFNAME);
345 rschregle 2.11 shmFile = mkstemp(shmFname);
346    
347 rschregle 2.13 if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0)
348     error(SYSTEM, "failed shared mem init in distribPhotons");
349 rschregle 2.11
350     photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE,
351     MAP_SHARED, shmFile, 0);
352    
353     if (photonCnt == MAP_FAILED)
354 rschregle 2.13 error(SYSTEM, "failed mapping shared memory in distribPhotons");
355     #else
356     /* Allocate photon counters statically on Windoze */
357     if (!(photonCnt = malloc(sizeof(PhotonCnt))))
358     error(SYSTEM, "failed trivial malloc in distribPhotons");
359     photonCnt -> numEmitted = photonCnt -> numComplete = 0;
360     #endif /* NIX */
361    
362     if (verbose) {
363     sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
364    
365     if (photonPorts) {
366     sprintf(errmsg2, " via %d ports", numPhotonPorts);
367     strcat(errmsg, errmsg2);
368     }
369    
370     strcat(errmsg, "\n");
371     eputs(errmsg);
372     }
373 greg 2.1
374     /* ===================================================================
375     * FLUX INTEGRATION - Get total photon flux from light sources
376     * =================================================================== */
377 rschregle 2.11 for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
378 greg 2.1 unsigned portCnt = 0;
379     emap.src = source + srcIdx;
380    
381 rschregle 2.11 do { /* Need at least one iteration if no ports! */
382 greg 2.1 emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
383     : NULL;
384     photonPartition [emap.src -> so -> otype] (&emap);
385    
386 rschregle 2.13 if (verbose) {
387     sprintf(errmsg, "\tIntegrating flux from source %s ",
388 greg 2.1 source [srcIdx].so -> oname);
389 rschregle 2.13
390 greg 2.1 if (emap.port) {
391     sprintf(errmsg2, "via port %s ",
392     photonPorts [portCnt].so -> oname);
393     strcat(errmsg, errmsg2);
394     }
395 rschregle 2.13
396     sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
397 greg 2.1 strcat(errmsg, errmsg2);
398     eputs(errmsg);
399 rschregle 2.13 #if NIX
400 greg 2.1 fflush(stderr);
401 rschregle 2.13 #endif
402 greg 2.1 }
403    
404     for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
405     emap.partitionCnt++) {
406     initPhotonEmission(&emap, pdfSamples);
407     totalFlux += colorAvg(emap.partFlux);
408     }
409    
410     portCnt++;
411     } while (portCnt < numPhotonPorts);
412     }
413    
414     if (totalFlux < FTINY)
415     error(USER, "zero flux from light sources");
416 rschregle 2.13
417     /* Record start time for progress reports */
418     repStartTime = time(NULL);
419    
420     if (verbose) {
421     sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
422     eputs(errmsg);
423     }
424 greg 2.1
425 rschregle 2.11 /* MAIN LOOP */
426     for (proc = 0; proc < numProc; proc++) {
427 rschregle 2.13 #if NIX
428 rschregle 2.11 if (!(pid = fork())) {
429 rschregle 2.13 /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */
430     #else
431     if (1) {
432     /* No subprocess under Windoze */
433     #endif
434     /* Local photon counters for this subprocess */
435 rschregle 2.11 unsigned passCnt = 0, prePassCnt = 0;
436     unsigned long lastNumPhotons [NUM_PMAP_TYPES];
437     unsigned long localNumEmitted = 0; /* Num photons emitted by this
438     subprocess alone */
439 greg 2.1
440 rschregle 2.11 /* Seed RNGs from PID for decorellated photon distribution */
441     pmapSeed(randSeed + proc, partState);
442 rschregle 2.13 pmapSeed(randSeed + (proc + 1) % numProc, emitState);
443     pmapSeed(randSeed + (proc + 2) % numProc, cntState);
444     pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
445     pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
446     pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
447 greg 2.1
448     for (t = 0; t < NUM_PMAP_TYPES; t++)
449 rschregle 2.11 lastNumPhotons [t] = 0;
450    
451     /* =============================================================
452     * 2-PASS PHOTON DISTRIBUTION
453     * Pass 1 (pre): emit fraction of target photon count
454     * Pass 2 (main): based on outcome of pass 1, estimate remaining
455     * number of photons to emit to approximate target
456     * count
457     * ============================================================= */
458 greg 2.1 do {
459 rschregle 2.11 double numEmit;
460 greg 2.1
461 rschregle 2.11 if (!passCnt) {
462     /* INIT PASS 1 */
463     /* Skip if no photons contributed after sufficient
464     * iterations; make it clear to user which photon maps are
465     * missing so (s)he can check geometry and materials */
466     if (++prePassCnt > maxPreDistrib) {
467 rschregle 2.13 sprintf(errmsg, "proc %d: too many prepasses", proc);
468 rschregle 2.11
469     for (t = 0; t < NUM_PMAP_TYPES; t++)
470     if (pmaps [t] && !pmaps [t] -> numPhotons) {
471     sprintf(errmsg2, ", no %s photons stored",
472     pmapName [t]);
473     strcat(errmsg, errmsg2);
474     }
475    
476     error(USER, errmsg);
477     break;
478 greg 2.1 }
479 rschregle 2.11
480     /* Num to emit is fraction of minimum target count */
481     numEmit = FHUGE;
482 greg 2.1
483 rschregle 2.11 for (t = 0; t < NUM_PMAP_TYPES; t++)
484     if (pmaps [t])
485     numEmit = min(pmaps [t] -> distribTarget, numEmit);
486    
487     numEmit *= preDistrib;
488 greg 2.1 }
489 rschregle 2.11 else {
490     /* INIT PASS 2 */
491     /* Based on the outcome of the predistribution we can now
492     * estimate how many more photons we have to emit for each
493     * photon map to meet its respective target count. This
494     * value is clamped to 0 in case the target has already been
495     * exceeded in the pass 1. */
496     double maxDistribRatio = 0;
497    
498     /* Set the distribution ratio for each map; this indicates
499     * how many photons of each respective type are stored per
500     * emitted photon, and is used as probability for storing a
501     * photon by newPhoton(). Since this biases the photon
502     * density, newPhoton() promotes the flux of stored photons
503     * to compensate. */
504     for (t = 0; t < NUM_PMAP_TYPES; t++)
505     if ((pm = pmaps [t])) {
506     pm -> distribRatio = (double)pm -> distribTarget /
507     pm -> numPhotons - 1;
508    
509     /* Check if photon map "overflowed", i.e. exceeded its
510     * target count in the prepass; correcting the photon
511     * flux via the distribution ratio is no longer
512     * possible, as no more photons of this type will be
513     * stored, so notify the user rather than deliver
514     * incorrect results. In future we should handle this
515     * more intelligently by using the photonFlux in each
516     * photon map to individually correct the flux after
517     * distribution. */
518     if (pm -> distribRatio <= FTINY) {
519     sprintf(errmsg, "%s photon map overflow in "
520     "prepass, reduce -apD", pmapName [t]);
521     error(INTERNAL, errmsg);
522     }
523    
524     maxDistribRatio = max(pm -> distribRatio,
525     maxDistribRatio);
526     }
527 greg 2.1
528 rschregle 2.11 /* Normalise distribution ratios and calculate number of
529     * photons to emit in main pass */
530     for (t = 0; t < NUM_PMAP_TYPES; t++)
531     if ((pm = pmaps [t]))
532     pm -> distribRatio /= maxDistribRatio;
533    
534     if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY)
535     /* No photons left to distribute in main pass */
536     break;
537     }
538    
539 rschregle 2.13 /* Update shared completion counter for progreport by parent */
540 rschregle 2.11 photonCnt -> numComplete += numEmit;
541    
542     /* PHOTON DISTRIBUTION LOOP */
543     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
544     unsigned portCnt = 0;
545     emap.src = source + srcIdx;
546    
547     do { /* Need at least one iteration if no ports! */
548     emap.port = emap.src -> sflags & SDISTANT
549     ? photonPorts + portCnt : NULL;
550     photonPartition [emap.src -> so -> otype] (&emap);
551    
552 rschregle 2.13 if (verbose && !proc) {
553     /* Output from subproc 0 only to avoid race condition
554     * on console I/O */
555 rschregle 2.11 if (!passCnt)
556 rschregle 2.13 sprintf(errmsg, "\tPREPASS %d on source %s ",
557 rschregle 2.11 prePassCnt, source [srcIdx].so -> oname);
558     else
559 rschregle 2.13 sprintf(errmsg, "\tMAIN PASS on source %s ",
560 rschregle 2.11 source [srcIdx].so -> oname);
561 rschregle 2.13
562 rschregle 2.11 if (emap.port) {
563     sprintf(errmsg2, "via port %s ",
564     photonPorts [portCnt].so -> oname);
565     strcat(errmsg, errmsg2);
566     }
567 rschregle 2.13
568 rschregle 2.11 sprintf(errmsg2, "(%lu partitions)\n",
569     emap.numPartitions);
570     strcat(errmsg, errmsg2);
571     eputs(errmsg);
572 rschregle 2.13 #if NIX
573 rschregle 2.11 fflush(stderr);
574 rschregle 2.13 #endif
575 rschregle 2.11 }
576 greg 2.1
577 rschregle 2.11 for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
578     emap.partitionCnt++) {
579     double partNumEmit;
580     unsigned long partEmitCnt;
581    
582     /* Get photon origin within current source partishunn
583     * and build emission map */
584     photonOrigin [emap.src -> so -> otype] (&emap);
585     initPhotonEmission(&emap, pdfSamples);
586    
587     /* Number of photons to emit from ziss partishunn --
588     * proportional to flux; photon ray weight and scalar
589 rschregle 2.13 * flux are uniform (latter only varying in RGB). */
590 rschregle 2.11 partNumEmit = numEmit * colorAvg(emap.partFlux) /
591     totalFlux;
592     partEmitCnt = (unsigned long)partNumEmit;
593    
594     /* Probabilistically account for fractional photons */
595     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
596     partEmitCnt++;
597    
598     /* Update local and shared (global) emission counter */
599     photonCnt -> numEmitted += partEmitCnt;
600     localNumEmitted += partEmitCnt;
601    
602     /* Integer counter avoids FP rounding errors during
603     * iteration */
604     while (partEmitCnt--) {
605     RAY photonRay;
606    
607     /* Emit photon based on PDF and trace through scene
608     * until absorbed/leaked */
609     emitPhoton(&emap, &photonRay);
610 rschregle 2.13 #if 1
611     if (emap.port)
612     /* !!! PHOTON PORT REJECTION SAMPLING HACK: set
613     * !!! photon port as fake hit object for
614     * !!! primary ray to check for intersection in
615     * !!! tracePhoton() */
616     photonRay.ro = emap.port -> so;
617     #endif
618 rschregle 2.11 tracePhoton(&photonRay);
619     }
620 rschregle 2.13
621 rschregle 2.11 /* Update shared global photon count for each pmap */
622     for (t = 0; t < NUM_PMAP_TYPES; t++)
623     if (pmaps [t]) {
624     photonCnt -> numPhotons [t] +=
625     pmaps [t] -> numPhotons - lastNumPhotons [t];
626     lastNumPhotons [t] = pmaps [t] -> numPhotons;
627     }
628 rschregle 2.13 #if !NIX
629     /* Synchronous progress report on Windoze */
630     if (!proc && photonRepTime > 0 &&
631     time(NULL) >= repLastTime + photonRepTime) {
632     repEmitted = repProgress = photonCnt -> numEmitted;
633     repComplete = photonCnt -> numComplete;
634     pmapDistribReport();
635     }
636     #endif
637 rschregle 2.11 }
638 greg 2.1
639 rschregle 2.11 portCnt++;
640     } while (portCnt < numPhotonPorts);
641     }
642    
643     for (t = 0; t < NUM_PMAP_TYPES; t++)
644     if (pmaps [t] && !pmaps [t] -> numPhotons) {
645     /* Double preDistrib in case a photon map is empty and
646     * redo pass 1 --> possibility of infinite loop for
647     * pathological scenes (e.g. absorbing materials) */
648     preDistrib *= 2;
649     break;
650 greg 2.1 }
651 rschregle 2.11
652 rschregle 2.13 if (t >= NUM_PMAP_TYPES)
653 rschregle 2.11 /* No empty photon maps found; now do pass 2 */
654     passCnt++;
655     } while (passCnt < 2);
656    
657 rschregle 2.13 /* Flush heap buffa for every pmap one final time;
658     * avoids potential data corruption! */
659 rschregle 2.11 for (t = 0; t < NUM_PMAP_TYPES; t++)
660     if (pmaps [t]) {
661     flushPhotonHeap(pmaps [t]);
662 rschregle 2.13 /* Heap file closed automatically on exit
663     fclose(pmaps [t] -> heap); */
664 rschregle 2.11 #ifdef DEBUG_PMAP
665 rschregle 2.13 sprintf(errmsg, "Proc %d: total %ld photons\n", proc,
666 rschregle 2.11 pmaps [t] -> numPhotons);
667     eputs(errmsg);
668     #endif
669     }
670 rschregle 2.13 #if NIX
671     /* Terminate subprocess */
672 rschregle 2.11 exit(0);
673 rschregle 2.13 #endif
674 greg 2.1 }
675 rschregle 2.11 else if (pid < 0)
676     error(SYSTEM, "failed to fork subprocess in distribPhotons");
677     }
678    
679 rschregle 2.13 #if NIX
680 rschregle 2.11 /* PARENT PROCESS CONTINUES HERE */
681     #ifdef SIGCONT
682 rschregle 2.13 /* Enable progress report signal handler */
683 rschregle 2.11 signal(SIGCONT, pmapDistribReport);
684 rschregle 2.13 #endif
685     /* Wait for subprocesses complete while reporting progress */
686 rschregle 2.11 proc = numProc;
687     while (proc) {
688     while (waitpid(-1, &stat, WNOHANG) > 0) {
689     /* Subprocess exited; check status */
690     if (!WIFEXITED(stat) || WEXITSTATUS(stat))
691     error(USER, "failed photon distribution");
692    
693     --proc;
694     }
695    
696     /* Nod off for a bit and update progress */
697     sleep(1);
698 rschregle 2.13
699     /* Asynchronous progress report from shared subprocess counters */
700 rschregle 2.11 repEmitted = repProgress = photonCnt -> numEmitted;
701 rschregle 2.13 repComplete = photonCnt -> numComplete;
702 rschregle 2.11
703 rschregle 2.13 repProgress = repComplete = 0;
704 greg 2.1 for (t = 0; t < NUM_PMAP_TYPES; t++)
705 rschregle 2.11 if ((pm = pmaps [t])) {
706     /* Get global photon count from shmem updated by subprocs */
707 rschregle 2.13 repProgress += pm -> numPhotons = photonCnt -> numPhotons [t];
708     repComplete += pm -> distribTarget;
709 greg 2.1 }
710 rschregle 2.13 repComplete *= numProc;
711 rschregle 2.11
712     if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
713     pmapDistribReport();
714     #ifdef SIGCONT
715     else signal(SIGCONT, pmapDistribReport);
716     #endif
717     }
718 rschregle 2.13 #endif /* NIX */
719 greg 2.1
720     /* ===================================================================
721 rschregle 2.11 * POST-DISTRIBUTION - Set photon flux and build data struct for photon
722     * storage, etc.
723 greg 2.1 * =================================================================== */
724 rschregle 2.11 #ifdef SIGCONT
725 rschregle 2.13 /* Reset signal handler */
726 rschregle 2.11 signal(SIGCONT, SIG_DFL);
727     #endif
728 greg 2.1 free(emap.samples);
729    
730 rschregle 2.13 /* Set photon flux */
731 rschregle 2.11 totalFlux /= photonCnt -> numEmitted;
732 rschregle 2.13 #if NIX
733 rschregle 2.11 /* Photon counters no longer needed, unmap shared memory */
734     munmap(photonCnt, sizeof(*photonCnt));
735     close(shmFile);
736 rschregle 2.13 unlink(shmFname);
737 rschregle 2.11 #else
738 rschregle 2.13 free(photonCnt);
739 rschregle 2.11 #endif
740 rschregle 2.13 if (verbose)
741     eputs("\n");
742    
743 greg 2.1 for (t = 0; t < NUM_PMAP_TYPES; t++)
744 rschregle 2.8 if (pmaps [t]) {
745 rschregle 2.13 if (verbose) {
746     sprintf(errmsg, "Building %s photon map\n", pmapName [t]);
747 greg 2.1 eputs(errmsg);
748 rschregle 2.13 #if NIX
749 greg 2.1 fflush(stderr);
750 rschregle 2.13 #endif
751 greg 2.1 }
752 rschregle 2.11
753     /* Build underlying data structure; heap is destroyed */
754     buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
755 greg 2.1 }
756 rschregle 2.13
757 greg 2.1 /* Precompute photon irradiance if necessary */
758 rschregle 2.13 if (preCompPmap) {
759     if (verbose)
760     eputs("\n");
761 greg 2.1 preComputeGlobal(preCompPmap);
762 rschregle 2.13 }
763    
764     if (verbose)
765     eputs("\n");
766 greg 2.1 }