Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmap.c
fab831a1d6fdpmapmerge

Last Change
Raw File

pmap.c
View Options

#ifndef lint
static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $";
#endif
/*
======================================================================
Photon map main module
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Fraunhofer Institute for Solar Energy Systems,
supported by the German Research Foundation
(DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS")
(c) Lucerne University of Applied Sciences and Arts,
supported by the Swiss National Science Foundation
(SNSF #147053, "Daylight Redirecting Components")
(c) Tokyo University of Science,
supported by the JSPS Grants-in-Aid for Scientific Research
(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
======================================================================
$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $
*/
#include "pmap.h"
#include "pmapmat.h"
#include "pmapsrc.h"
#include "pmaprand.h"
#include "pmaproi.h"
#include "pmapio.h"
#include "pmapdens.h"
#include "pmapbias.h"
#include "pmapdiag.h"
#include "pmutil.h"
#include "otypes.h"
#include "otspecial.h"
#include <time.h>
#if NIX
#include <sys/stat.h>
#include <sys/mman.h>
#include <sys/wait.h>
#include <sys/file.h>
#endif
static int photonParticipate (RAY *ray)
/* Trace photon through participating medium. Returns 1 if passed through,
or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */
{
int i;
RREAL xi1, cosTheta, phi, du, dv, rmax0;
const RREAL cext = colorAvg(ray -> cext), invCext = 1 / cext,
albedo = colorAvg(ray -> albedo), invAlbedo = 1 / albedo,
gecc = ray -> gecc, gecc2 = sqr(gecc);
FVECT u, v;
COLOR cvext;
/* Save incoming ray's max distance; this is required since rmax is
set to the mean free distance for each path segment (see below) */
rmax0 = ray -> rmax;
/* Mean free distance until interaction with medium */
ray -> rmax = -log(pmapRandom(mediumState)) / cext;
while (!localhit(ray, &thescene)) {
/* No geometry hit within mean free dist ray->rmax; photon interacts
with medium */
if (!incube(&thescene, ray -> rop)) {
/* Terminate photon if it has leaked from the scene */
#ifdef DEBUG_PMAP_LEAKED
fprintf(
stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n",
ray -> rop [0], ray -> rop [1], ray -> rop [2]
);
#endif
return 0;
}
/* RGB attenuation over mean free distance */
setcolor(cvext,
exp(-ray -> rmax * ray -> cext [0]),
exp(-ray -> rmax * ray -> cext [1]),
exp(-ray -> rmax * ray -> cext [2])
);
/* Modify ray color and normalise */
multcolor(ray -> rcol, cvext);
colorNorm(ray -> rcol);
/* Set ray origin for next interaction */
VCOPY(ray -> rorg, ray -> rop);
/* Update ray's path length; terminate photon if maximum reached.
See also photonRay(). */
rmax0 -= ray -> rot;
if (photonMaxPathDist > 0 && rmax0 < 0)
return 0;
/* Store volume photon if nonzero albedo; this also accounts for
direct inscattering from light sources */
if (albedo > FTINY) {
/* Set ray's path length to store with photon */
ray -> rmax = rmax0;
#ifdef PMAP_PHOTONFLOW
if (lightFlowPhotonMapping) {
/* Temporarily correct normalised flux for lightflow photons
based on extinction (= mean num photons deposited per unit
path length). */
scalecolor(ray -> rcol, invCext);
newPhoton(lightFlowPmap, ray);
/* Undo correction for next iteration after storing photon */
scalecolor(ray -> rcol, cext);
}
else
#endif
newPhoton(volumePmap, ray);
}
/* Colour bleeding without attenuation (?) */
multcolor(ray -> rcol, ray -> albedo);
scalecolor(ray -> rcol, invAlbedo);
if (pmapRandom(rouletteState) < albedo) {
/* Scatter photon */
xi1 = pmapRandom(scatterState);
cosTheta = ray -> gecc <= FTINY
? 2 * xi1 - 1
: 0.5 / gecc * (
1 + gecc2 - sqr((1 - gecc2) / (1 + gecc * (2 * xi1 - 1)))
);
phi = 2 * PI * pmapRandom(scatterState);
du = dv = sqrt(1 - sqr(cosTheta)); /* sin(theta) */
du *= cos(phi);
dv *= sin(phi);
/* Get axes u & v perpendicular to photon direction */
i = 0;
do {
v [0] = v [1] = v [2] = 0;
v [i++] = 1;
fcross(u, v, ray -> rdir);
} while (normalize(u) < FTINY);
fcross(v, ray -> rdir, u);
for (i = 0; i < 3; i++)
ray -> rdir [i] = du * u [i] + dv * v [i] +
cosTheta * ray -> rdir [i];
}
ray -> rlvl++;
ray -> rmax = -log(pmapRandom(mediumState)) / cext;
/* fprintf(stderr, "%.3f\n", ray -> rmax); */
}
/* Passed through medium until intersecting local object */
setcolor(cvext,
exp(-ray -> rot * ray -> cext [0]),
exp(-ray -> rot * ray -> cext [1]),
exp(-ray -> rot * ray -> cext [2])
);
/* Modify ray color and normalise */
multcolor(ray -> rcol, cvext);
colorNorm(ray -> rcol);
/* Restore outgoing ray's max distance */
ray -> rmax = rmax0;
return 1;
}
void tracePhoton (RAY *ray)
/* Follow photon as it bounces around... */
{
long mod;
OBJREC *mat, *port = NULL;
if (!ray -> parent) {
/* !!! PHOTON PORT REJECTION SAMPLING HACK: get photon port for
* !!! primary ray from ray -> ro, then reset the latter to NULL so
* !!! as not to interfere with localhit() */
port = ray -> ro;
ray -> ro = NULL;
}
if (ray -> rlvl > photonMaxBounce) {
#ifdef PMAP_RUNAWAY_WARN
error(WARNING, "runaway photon!");
#endif
return;
}
if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
return;
/* NOTE: localhit() limits intersections to the aft plane's distance
ray->rmax. This mechanism is (ab)used here to limit the photon path
length by initialising ray->rmax to photonMaxPathDist; see emitPhoton().
With unlimited path length (photonMaxPathDist=0), ray->rmax becomes
negative (see photonRay()), and the aft plane has no effect. */
if (localhit(ray, &thescene)) {
mod = ray -> ro -> omod;
if (port && ray -> ro != port) {
/* Terminate photon if emitted from port without intersecting it */
#ifdef PMAP_PORTREJECT_WARN
sprintf(errmsg, "photon outside port %s", ray -> ro -> oname);
error(WARNING, errmsg);
#endif
return;
}
if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) {
/* Transfer ray if modifier is VOID or clipped within antimatta */
RAY tray;
photonRay(ray, &tray, PMAP_XFER, NULL);
tracePhoton(&tray);
}
else {
/* Scatter for modifier material */
mat = objptr(mod);
photonScatter [mat -> otype] (mat, ray);
}
}
}
static void preComputeGlobal (PhotonMap *pmap)
/* Precompute irradiance from global photons for final gathering for
a random subset of finalGather * pmap -> numPhotons photons, and builds
the photon map, discarding the original photons. */
/* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */
{
unsigned long i, numPreComp;
unsigned j;
PhotonIdx pIdx;
Photon photon;
RAY ray;
PhotonMap nuPmap;
repComplete = numPreComp = finalGather * pmap -> numPhotons;
if (verbose) {
sprintf(errmsg,
"\nPrecomputing irradiance for %ld global photons\n", numPreComp
);
eputs(errmsg);
#if NIX
fflush(stderr);
#endif
}
/* Copy photon map for precomputed photons */
memcpy(&nuPmap, pmap, sizeof(PhotonMap));
/* Zero counters, init new heap and extents */
nuPmap.numPhotons = 0;
initPhotonHeap(&nuPmap);
for (j = 0; j < 3; j++) {
nuPmap.minPos [j] = FHUGE;
nuPmap.maxPos [j] = -FHUGE;
}
/* Record start time, baby */
repStartTime = time(NULL);
#ifdef SIGCONT
signal(SIGCONT, pmapPreCompReport);
#endif
repProgress = 0;
photonRay(NULL, &ray, PRIMARY, NULL);
ray.ro = NULL;
for (i = 0; i < numPreComp; i++) {
/* Get random photon from stratified distribution in source heap to
* avoid duplicates and clustering */
pIdx = firstPhoton(pmap) +
(unsigned long)((i + pmapRandom(pmap -> randState)) /
finalGather);
getPhoton(pmap, pIdx, &photon);
/* Init dummy photon ray with intersection at photon position */
VCOPY(ray.rop, photon.pos);
for (j = 0; j < 3; j++)
ray.ron [j] = photon.norm [j] / 127.0;
/* Get density estimate at photon position */
photonDensity(pmap, &ray, ray.rcol);
/* Append photon to new heap from ray */
newPhoton(&nuPmap, &ray);
/* Update progress */
repProgress++;
if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
pmapPreCompReport();
#ifdef SIGCONT
else signal(SIGCONT, pmapPreCompReport);
#endif
}
/* Flush heap */
flushPhotonHeap(&nuPmap);
#ifdef SIGCONT
signal(SIGCONT, SIG_DFL);
#endif
/* Trash original pmap, replace with precomputed one */
deletePhotons(pmap);
memcpy(pmap, &nuPmap, sizeof(PhotonMap));
if (verbose) {
eputs("\nRebuilding precomputed photon map\n");
#if NIX
fflush(stderr);
#endif
}
/* Rebuild underlying data structure, destroying heap */
buildPhotonMap(pmap, NULL, NULL, 1);
}
typedef struct {
unsigned long numPhotons [NUM_PMAP_TYPES],
numEmitted, numComplete;
} PhotonCnt;
void distribPhotons (PhotonMap **pmaps, unsigned numProc)
{
EmissionMap emap;
char errmsg2 [128], shmFname [PMAP_TMPFNLEN];
unsigned t, srcIdx, proc;
double totalFlux = 0;
int shmFile, stat, pid;
PhotonMap *pm;
PhotonCnt *photonCnt;
pid_t procPids [PMAP_MAXPROC];
for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++);
if (t >= NUM_PMAP_TYPES)
error(USER, "no photon maps defined in distribPhotons");
if (!nsources)
error(USER, "no light sources in distribPhotons");
/* ===================================================================
* INITIALISATION - Set up emission and scattering funcs
* =================================================================== */
emap.samples = NULL;
emap.maxPartitions = MAXSPART;
emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
if (!emap.partitions)
error(INTERNAL, "can't allocate source partitions in distribPhotons");
/* Initialise all defined photon maps */
for (t = 0; t < NUM_PMAP_TYPES; t++)
if (pmaps [t]) {
initPhotonMap(pmaps [t], t);
/* Open photon heapfile */
initPhotonHeap(pmaps [t]);
/* Per-subprocess target count */
pmaps [t] -> distribTarget /= numProc;
if (!pmaps [t] -> distribTarget)
error(INTERNAL, "no photons to distribute in distribPhotons");
}
initPhotonEmissionFuncs();
initPhotonScatterFuncs();
/* Get photon ports from modifier list */
getPhotonPorts(photonPortList);
/* Get photon sensor modifiers */
getPhotonSensors(photonSensorList);
/* Get polyhedral regions of interest */
getPolyROIs(pmapROImodList);
#if NIX
/* Set up shared mem for photon counters (zeroed by ftruncate) */
strcpy(shmFname, PMAP_TMPFNAME);
shmFile = mkstemp(shmFname);
if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0)
error(SYSTEM, "failed shared mem init in distribPhotons");
photonCnt = mmap(NULL, sizeof(*photonCnt),
PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0
);
if (photonCnt == MAP_FAILED)
error(SYSTEM, "failed mapping shared memory in distribPhotons");
#else
/* Allocate photon counters statically on Windoze */
if (!(photonCnt = malloc(sizeof(PhotonCnt))))
error(SYSTEM, "failed trivial malloc in distribPhotons");
photonCnt -> numEmitted = photonCnt -> numComplete = 0;
#endif /* NIX */
if (verbose) {
sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
if (photonPorts) {
sprintf(errmsg2, " via %d ports", numPhotonPorts);
strcat(errmsg, errmsg2);
}
strcat(errmsg, "\n");
eputs(errmsg);
}
/* ===================================================================
* FLUX INTEGRATION - Get total photon flux from light sources
* =================================================================== */
for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
unsigned portCnt = 0;
emap.src = source + srcIdx;
do {
/* Set next photon port if defined; note this loop iterates once if
* no ports are defined. */
emap.port = emap.src -> sflags & SDISTANT
? photonPorts + portCnt
: NULL;
photonPartition [emap.src -> so -> otype] (&emap);
if (verbose) {
sprintf(errmsg, "\tIntegrating flux from source %s ",
source [srcIdx].so -> oname
);
if (emap.port) {
sprintf(errmsg2, "via port %s ",
photonPorts [portCnt].so -> oname
);
strcat(errmsg, errmsg2);
}
sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
strcat(errmsg, errmsg2);
eputs(errmsg);
#if NIX
fflush(stderr);
#endif
}
for (emap.partitionCnt = 0;
emap.partitionCnt < emap.numPartitions;
emap.partitionCnt++
) {
initPhotonEmission(&emap, pdfSamples);
totalFlux += colorAvg(emap.partFlux);
}
portCnt++;
} while (portCnt < numPhotonPorts);
}
if (totalFlux < FTINY)
error(USER, "zero flux from light sources");
/* Record start time for progress reports */
repStartTime = time(NULL);
if (verbose) {
sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
eputs(errmsg);
}
/* MAIN LOOP */
for (proc = 0; proc < numProc; proc++) {
#if NIX
if (!(pid = fork())) {
/* SUBPROCESS ENTERS HERE; open and mmapped files inherited */
#else
if (1) {
/* No subprocess under Windoze */
#endif
/* Local photon counters for this subprocess */
unsigned passCnt = 0, prePassCnt = 0;
unsigned long preIncPhotonCnt, lastNumPhotons [NUM_PMAP_TYPES],
numEmitted = 0, lastNumEmitted = 0;
/* Seed RNGs from PID for decorellated photon distribution */
pmapSeed(randSeed + proc, partState);
pmapSeed(randSeed + (proc + 1) % numProc, emitState);
pmapSeed(randSeed + (proc + 2) % numProc, cntState);
pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
#ifdef DEBUG_PMAP
/* Output child process PID after random delay to prevent corrupted
* console output due to race condition */
usleep(1e6 * pmapRandom(rouletteState));
fprintf(stderr,
"Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n",
proc, getpid()
);
/* Allow time for debugger to attach to child process */
sleep(10);
#endif
for (t = 0; t < NUM_PMAP_TYPES; t++)
lastNumPhotons [t] = 0;
/* =============================================================
* 2-PASS PHOTON DISTRIBUTION
* Pass 1 (pre): emit fraction of target photon count
* Pass 2 (main): based on outcome of pass 1, estimate remaining
* number of photons to emit to approximate target
* count
* ============================================================= */
do {
double numEmit;
if (!passCnt) {
/* INIT PASS 1 */
/* Skip if no photons contributed after sufficient
* iterations; make it clear to user which photon maps are
* missing so (s)he can check geometry and materials */
if (++prePassCnt > maxPreDistrib) {
sprintf(errmsg, "proc %d: too many prepasses", proc);
for (t = 0; t < NUM_PMAP_TYPES; t++)
if (pmaps [t] && !pmaps [t] -> numPhotons) {
sprintf(errmsg2, ", no %s photons stored",
pmapName [t]
);
strcat(errmsg, errmsg2);
}
sprintf(errmsg2, "; increase -apD");
strcat(errmsg, errmsg2);
error(USER, errmsg);
break;
}
/* Num to emit is fraction of minimum target count */
numEmit = FHUGE;
for (t = 0; t < NUM_PMAP_TYPES; t++)
if (pmaps [t])
numEmit = min(pmaps [t] -> distribTarget, numEmit);
numEmit *= preDistrib;
}
else {
/* INIT PASS 2 */
/* Based on the outcome of the predistribution we can now
* estimate how many more photons we have to emit for each
* photon map to meet its respective target count. This
* value is clamped to 0 in case the target has already been
* exceeded in the pass 1. */
double maxDistribRatio = 0;
/* Set the distribution ratio for each map; this indicates
* how many photons of each respective type are stored per
* emitted photon, and is used as probability for storing a
* photon by newPhoton(). Since this biases the photon
* density, newPhoton() promotes the flux of stored photons
* to compensate. */
for (t = 0; t < NUM_PMAP_TYPES; t++)
if ((pm = pmaps [t])) {
pm -> distribRatio =
(double)pm -> distribTarget / pm -> numPhotons - 1;
/* Check if photon map "overflowed", i.e. exceeded its
* target count in the prepass; correcting the photon
* flux via the distribution ratio is no longer
* possible, as no more photons of this type will be
* stored, so notify the user rather than deliver
* incorrect results. In future we should handle this
* more intelligently by using the photonFlux in each
* photon map to individually correct the flux after
* distribution. */
if (pm -> distribRatio <= FTINY) {
sprintf(errmsg,
"%s photon map overflow in prepass, reduce -apD ",
pmapName [t]
);
error(INTERNAL, errmsg);
}
maxDistribRatio = max(pm -> distribRatio,
maxDistribRatio
);
}
/* Normalise distribution ratios and calculate number of
* photons to emit in main pass */
for (t = 0; t < NUM_PMAP_TYPES; t++)
if ((pm = pmaps [t]))
pm -> distribRatio /= maxDistribRatio;
if ((numEmit = numEmitted * maxDistribRatio) < FTINY)
/* No photons left to distribute in main pass */
break;
}
/* PHOTON DISTRIBUTION LOOP */
for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
unsigned portCnt = 0;
emap.src = source + srcIdx;
do {
/* Set next photon port if defined; note this loop iterates
* once if no ports are defined. */
emap.port = emap.src -> sflags & SDISTANT
? photonPorts + portCnt
: NULL;
photonPartition [emap.src -> so -> otype] (&emap);
if (verbose && !proc) {
/* Output from subproc 0 only to avoid race condition
* on console I/O */
if (!passCnt)
sprintf(errmsg, "\tPREPASS %d on source %s ",
prePassCnt, source [srcIdx].so -> oname
);
else
sprintf(errmsg, "\tMAIN PASS on source %s ",
source [srcIdx].so -> oname
);
if (emap.port) {
sprintf(errmsg2, "via port %s ",
photonPorts [portCnt].so -> oname
);
strcat(errmsg, errmsg2);
}
sprintf(errmsg2, "(%lu partitions)\n",
emap.numPartitions
);
strcat(errmsg, errmsg2);
eputs(errmsg);
#if NIX
fflush(stderr);
#endif
}
for (emap.partitionCnt = 0;
emap.partitionCnt < emap.numPartitions;
emap.partitionCnt++
) {
double partNumEmit;
unsigned long partEmitCnt;
/* Get photon origin within current source partishunn
* and build emission map */
photonOrigin [emap.src -> so -> otype] (&emap);
initPhotonEmission(&emap, pdfSamples);
/* Number of photons to emit from ziss partishunn --
* proportional to flux; photon ray weight and scalar
* flux are uniform (latter only varying in RGB). */
partNumEmit = numEmit * colorAvg(emap.partFlux) /
totalFlux;
partEmitCnt = (unsigned long)partNumEmit;
/* Probabilistically account for fractional photons */
if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
partEmitCnt++;
/* Integer counter avoids FP rounding errors during
* iteration */
while (partEmitCnt--) {
RAY photonRay;
/* Emit photon based on PDF, skip if rejected. */
/* NOTE: rejection sampling skips incrementing the
* emission counter (see below), thus compensating
* for the rejected photons by increasing the photon
* flux in proportion to the lower effective
* emission count.
* BUG: THIS INTERFERES WITH THE PROGRESS COUNTER
* REPORTED TO THE PARENT, AND WITH THE
* PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL
* FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED!
*/
if (!emitPhoton(&emap, &photonRay))
continue;
/* Set subprocess ID in photonRay; extends path ID.
Used by newPhoton() as photon attributes. */
PMAP_SETRAYPROC(&photonRay, proc);
tracePhoton(&photonRay);
/* Update local emission counter*/
numEmitted++;
if (!(partEmitCnt && (numEmitted & PMAP_CNTUPDATE))) {
/* Update global counters shared with siblings in
intervals to reduce overhead & contention;
ALSO WHEN EMITTING FINAL PHOTON FROM THIS
PARTITION, OTHERWISE COUNTERS WILL SKIP LAST
UPDATE, BIASING THE PHOTON FLUX! */
/* Differentially increment number emitted using
* local counter numEmitted */
preIncPhotonCnt = photonCnt -> numEmitted;
photonCnt -> numEmitted +=
numEmitted - lastNumEmitted;
photonCnt -> numComplete +=
numEmitted - lastNumEmitted;
lastNumEmitted = numEmitted;
/* Check overflow using pre-increment values */
if (photonCnt -> numEmitted < preIncPhotonCnt) {
sprintf(errmsg, "photon emission counter "
"overflow (was: %ld, is: %ld)",
preIncPhotonCnt, photonCnt -> numEmitted
);
error(INTERNAL, errmsg);
}
/* Update global photon count for each pmap */
for (t = 0; t < NUM_PMAP_TYPES; t++) {
if (pmaps [t]) {
preIncPhotonCnt = photonCnt->numPhotons [t];
/* Differentially increment using local
counter lastNumPhotons */
photonCnt -> numPhotons [t] +=
pmaps [t] -> numPhotons -
lastNumPhotons [t];
lastNumPhotons [t] =
pmaps [t] -> numPhotons;
/* Check for photon counter overflow (this
* could only happen before an emission
* counter overflow if the scene has an
* absurdly high albedo and/or very dense
* geometry) */
if (photonCnt -> numPhotons [t] <
preIncPhotonCnt
) {
sprintf(errmsg, "photon counter "
"overflow (was: %ld, is: %ld)",
preIncPhotonCnt,
photonCnt -> numPhotons [t]
);
error(INTERNAL, errmsg);
}
}
}
}
}
#if !NIX
/* Synchronous progress report on Windoze */
if (!proc && photonRepTime > 0 &&
time(NULL) >= repLastTime + photonRepTime) {
repEmitted = repProgress = photonCnt -> numEmitted;
repComplete = photonCnt -> numComplete;
pmapDistribReport();
}
#endif
}
portCnt++;
} while (portCnt < numPhotonPorts);
}
for (t = 0; t < NUM_PMAP_TYPES; t++)
if (pmaps [t] && !pmaps [t] -> numPhotons) {
/* Double preDistrib in case a photon map is empty and
* redo pass 1 --> possibility of infinite loop for
* pathological scenes (e.g. absorbing materials) */
preDistrib *= 2;
break;
}
if (t >= NUM_PMAP_TYPES)
/* No empty photon maps found; now do pass 2 */
passCnt++;
} while (passCnt < 2);
/* Flush heap buffa for every pmap one final time;
* avoids potential data corruption! */
for (t = 0; t < NUM_PMAP_TYPES; t++)
if (pmaps [t]) {
flushPhotonHeap(pmaps [t]);
/* Heap file closed automatically on exit
fclose(pmaps [t] -> heap); */
#ifdef DEBUG_PMAP
sprintf(errmsg, "Proc %d: total %ld photons\n",
proc, pmaps [t] -> numPhotons
);
eputs(errmsg);
#endif
}
#if NIX
/* Terminate subprocess */
exit(0);
#endif
}
else {
/* PARENT PROCESS CONTINUES LOOP ITERATION HERE */
if (pid < 0)
error(SYSTEM, "failed to fork subprocess in distribPhotons");
else
/* Save child process IDs */
procPids [proc] = pid;
}
}
#if NIX
/* PARENT PROCESS CONTINUES HERE */
#ifdef SIGCONT
/* Enable progress report signal handler */
signal(SIGCONT, pmapDistribReport);
#endif
/* Wait for subprocesses to complete while reporting progress */
proc = numProc;
while (proc) {
while (waitpid(-1, &stat, WNOHANG) > 0) {
/* Subprocess exited; check status */
if (!WIFEXITED(stat) || WEXITSTATUS(stat)) {
/* Exited with error; terminate siblings, clean up & bail out */
for (proc = 0; proc < numProc; proc++)
kill(procPids [proc], SIGKILL);
/* Unmap shared memory, squish mapped file */
munmap(photonCnt, sizeof(*photonCnt));
close(shmFile);
unlink(shmFname);
error(USER, "failed photon distribution");
}
--proc;
}
/* Nod off for a bit and update progress */
sleep(1);
/* Asynchronous progress report from shared subprocess counters */
repEmitted = repProgress = photonCnt -> numEmitted;
repComplete = photonCnt -> numComplete;
repProgress = repComplete = 0;
for (t = 0; t < NUM_PMAP_TYPES; t++)
if ((pm = pmaps [t])) {
/* Get global photon count from shmem updated by subprocs */
repProgress += pm -> numPhotons = photonCnt -> numPhotons [t];
repComplete += pm -> distribTarget;
}
repComplete *= numProc;
if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
pmapDistribReport();
#ifdef SIGCONT
else signal(SIGCONT, pmapDistribReport);
#endif
}
#endif /* NIX */
/* ===================================================================
* POST-DISTRIBUTION - Set photon flux and build data struct for photon
* storage, etc.
* =================================================================== */
#ifdef SIGCONT
/* Reset signal handler */
signal(SIGCONT, SIG_DFL);
#endif
free(emap.samples);
printf("numEmitted/CNTUPDATE: %f\n",
photonCnt -> numEmitted / (float)PMAP_CNTUPDATE
);
/* Set photon flux */
totalFlux /= photonCnt -> numEmitted;
#if NIX
/* Photon counters no longer needed, unmap shared memory */
munmap(photonCnt, sizeof(*photonCnt));
close(shmFile);
unlink(shmFname);
#else
free(photonCnt);
#endif
if (verbose)
eputs("\n");
for (t = 0; t < NUM_PMAP_TYPES; t++)
if (pmaps [t]) {
if (verbose) {
sprintf(errmsg, "Building %s photon map\n", pmapName [t]);
eputs(errmsg);
#if NIX
fflush(stderr);
#endif
}
/* Build underlying data structure; heap is destroyed */
buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
}
/* Precompute photon irradiance if necessary */
if (preCompPmap) {
if (verbose)
eputs("\n");
preComputeGlobal(preCompPmap);
}
if (verbose)
eputs("\n");
}
c4science · Help