Homec4science
Diffusion RADIANCE Photon Map (fixSpecular) Browse

RadiancePmap/pmapcontrib.c
bd570b55223efixSpecular

Last Change
Raw File

pmapcontrib.c
View Options

#ifndef lint
static const char RCSid[] = "$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $";
#endif
/*
=========================================================================
Photon map routines for precomputed light source contributions; these
routines interface to mkpmap.
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Lucerne University of Applied Sciences and Arts,
supported by the Swiss National Science Foundation
(SNSF #147053, "Daylight Redirecting Components",
SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment")
=========================================================================
$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $
*/
#include "pmapcontrib.h"
#include "pmapmat.h"
#include "pmapsrc.h"
#include "pmaprand.h"
#include "pmaproi.h"
#include "pmapio.h"
#include "pmapdiag.h"
#include "wavelet2.h"
#include "otypes.h"
#include "otspecial.h"
#if NIX
#include <sys/mman.h>
#include <sys/wait.h>
#endif
extern void SDdisk2square(double sq[2], double diskx, double disky);
MODCONT *addContribModifier (
LUTAB *contribTab, unsigned *numContribs,
char *mod, char *binParm, char *binVal, int binCnt
)
/* Add light source modifier mod to contribution lookup table contribTab,
and update numContribs. Return initialised contribution data for this
modifier. */
{
LUENT *lutEntry = lu_find(contribTab, mod);
MODCONT *contrib;
EPNODE *eBinVal;
if (lutEntry -> data) {
/* Reject duplicate modifiers */
sprintf(errmsg, "duplicate light source modifier %s", mod);
error(USER, errmsg);
}
if (*numContribs >= MAXMODLIST) {
sprintf(errmsg, "too many modifiers (max. %d)", MAXMODLIST);
error(INTERNAL, errmsg);
}
lutEntry -> key = mod;
if (!binVal)
/* Fall back to single bin if no value specified */
binVal = "0";
eBinVal = eparse(binVal);
if (eBinVal -> type == NUM) {
/* Bin value is constant */
binCnt = (int)(evalue(eBinVal) + 1.5);
if (binCnt != 1) {
sprintf(errmsg, "invalid non-zero constant for bin %s", binVal);
error(USER, errmsg);
}
}
else if (binCnt <= 0) {
sprintf(errmsg, "undefined/invalid bin count for modifier %s", mod);
error(USER, errmsg);
}
/* Allocate and init contributions */
contrib = (MODCONT *)malloc(
sizeof(MODCONT) + sizeof(DCOLOR) * (binCnt - 1)
);
if (!contrib)
error(SYSTEM, "out of memory in addContribModifier()");
contrib -> outspec = NULL;
contrib -> modname = mod;
contrib -> params = binParm;
contrib -> nbins = binCnt;
contrib -> binv = eBinVal;
contrib -> bin0 = 0;
memset(contrib -> cbin, 0, sizeof(DCOLOR) * binCnt);
lutEntry -> data = lutEntry -> data = (char *)contrib;
++(*numContribs);
return(contrib);
}
void addContribModfile (
LUTAB *contribTab, unsigned *numContribs,
char *modFile, char *binParm, char *binVal, int binCnt
)
/* Add light source modifiers from file modFile to contribution lookup table
* contribTab, and update numContribs.
* NOTE: This code is adapted from rcontrib.c */
{
char *mods [MAXMODLIST];
int i;
/* Find file and store strings */
i = wordfile(mods, MAXMODLIST, getpath(modFile, getrlibpath(), R_OK));
if (i < 0) {
sprintf(errmsg, "can't open modifier file %s", modFile);
error(SYSTEM, errmsg);
}
if (*numContribs + i >= MAXMODLIST - 1) {
sprintf(errmsg, "too many modifiers (max. %d) in file %s",
MAXMODLIST - 1, modFile
);
error(INTERNAL, errmsg);
}
for (i = 0; mods [i]; i++)
/* Add each modifier */
addContribModifier(contribTab, numContribs, mods [i],
binParm, binVal, binCnt
);
}
static unsigned logDim (unsigned size)
/* Return log2(sqrt(size)) = l, where size = (2^l)(2^l).
Is there a faster (and portable!) way? */
{
unsigned i, dim, sz;
if (!size)
return 0;
for (i = 0, dim = 1, sz = size; sz >>= 1; dim++);
return dim >> 1;
}
static int serialContribSourceBin (unsigned l, int x, int y)
/* Serialise 2D coordinates to 1D for array of size (2^l) x (2^l).
Returns -1 if coordinates invalid */
{
return x < 0 || y < 0
? -1
: (x << l) + y;
}
static void deserialContribSourceBin (unsigned l, int bin, int *x, int *y)
/* Deserialise 1D bin to 2D coordinates for array of size (2^l) x (2^l).
Returns -1 if bin invalid */
{
/* Obviously this is faster than integer division/modulo */
*x = bin < 0 ? -1 : bin >> l;
*y = bin < 0 ? -1 : bin & (1 << l) - 1;
}
static int contribSourceBin (LUTAB *contribs, const RAY *ray)
/* Map contribution source ray to its bin for light source ray -> rsrc,
using Shirley-Chiu disk-to-square mapping.
Return invalid bin -1 if mapping failed. */
{
const SRCREC *src;
const OBJREC *srcMod;
const MODCONT *srcCont;
int i, bin, scRHS;
FVECT scNorm, scUp;
RAY srcRay;
double diskx, disky, scCoord [2];
unsigned l, scDim, scBin [2];
/* Check we have a valid ray and contribution LUT */
if (!ray || !contribs)
return -1;
src = &source [ray -> rsrc];
srcMod = findmaterial(src -> so);
srcCont = (MODCONT*)lu_find(contribs, srcMod -> oname) -> data;
if (!srcCont)
/* We're not interested in this source (modifier not in contrib LUT) */
return -1;
/* Set up shadow ray pointing to source for disk2square mapping */
rayorigin(&srcRay, SHADOW, NULL, NULL);
srcRay.rsrc = ray -> rsrc;
VCOPY(srcRay.rorg, ray -> rop);
for (i = 0; i < 3; i++)
srcRay.rdir [i] = -ray -> rdir [i];
if (!(src -> sflags & SDISTANT
? sourcehit(&srcRay)
: (*ofun[srcMod -> otype].funp)(srcMod, &srcRay)
))
/* (Redundant?) sanity check for valid source ray? */
return -1;
worldfunc(RCCONTEXT, &srcRay);
set_eparams((char*)srcCont -> params);
/* Get Shirley-Chiu mapping orientation from function parameters */
scRHS = varvalue(PMAP_CONTRIB_SCRHS);
scNorm [0] = varvalue(PMAP_CONTRIB_SCNORMX);
scNorm [1] = varvalue(PMAP_CONTRIB_SCNORMY);
scNorm [2] = varvalue(PMAP_CONTRIB_SCNORMZ);
scUp [0] = varvalue(PMAP_CONTRIB_SCUPX);
scUp [1] = varvalue(PMAP_CONTRIB_SCUPY);
scUp [2] = varvalue(PMAP_CONTRIB_SCUPZ);
/* Apply Shirley-Chiu mapping of disk to square */
SDdisk2square(scCoord, diskx, disky);
l = logDim(srcCont -> nbins);
scDim = 1 << l;
/* Map Shirley-Chiu square coords to 2D bins and serialise */
scBin [0] = scCoord [0] * scDim + 0.5;
scBin [1] = scCoord [1] * scDim + 0.5;
if ((bin = serialContribSourceBin(l, scBin [0], scBin [1])) < 0) {
error(WARNING, "Ignoring invalid bin in contribSourceBin()");
return -1;
}
return bin;
}
static PhotonContribSourceIdx newPhotonContribSource (
PhotonMap *pmap, const RAY *contribSrcRay, FILE *contribSrcHeap
)
/* Add contribution source ray for emitted contribution photon and save
* light source index and binned direction. The current contribution source
* is stored in pmap -> lastContribSrc. If the previous contribution source
* spawned photons (i.e. has srcIdx >= 0), it's appended to contribSrcHeap.
* If contribSrcRay == NULL, the current contribution source is still
* flushed, but no new source is set. Returns updated contribution source
* counter pmap -> numContribSrc. */
{
if (!pmap || !contribSrcHeap)
return 0;
/* Check if last contribution source has spawned photons (srcIdx >= 0,
* see newPhoton()), in which case we save it to the heap file before
* clobbering it. (Note this is short-term storage, so we doan' need
* da portable I/O stuff here). */
if (pmap -> lastContribSrc.srcIdx >= 0) {
if (!fwrite(&pmap -> lastContribSrc,
sizeof(PhotonContribSource), 1, contribSrcHeap
))
error(SYSTEM,
"failed writing photon contrib source in newPhotonContribSource()"
);
pmap -> numContribSrc++;
if (!pmap -> numContribSrc ||
pmap -> numContribSrc > PMAP_MAXCONTRIBSRC
);
}
/* Mark this contribution source unused with a negative source index
until its path spawns a photon (see newPhoton()) */
pmap -> lastContribSrc.srcIdx = -1;
/* Map ray to bin in anticipation that this contribution source will
be used, since the ray will be lost once a photon is spawned */
pmap -> lastContribSrc.srcBin = contribSourceBin(
pmap -> contribTab, contribSrcRay
);
return pmap -> numContribSrc;
}
static PhotonContribSourceIdx buildContribSources (PhotonMap *pmap,
FILE **contribSrcHeap, char **contribSrcHeapFname,
PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps
)
/* Consolidate per-subprocess contribution sources heaps into array
* pmap -> contribSrc. Returns offset for contribution source index
* linearisation in pmap -> numContribSrc. The heap files in contribSrcHeap
* are closed on return. */
{
PhotonContribSourceIdx heapLen;
unsigned heap;
if (!pmap || !contribSrcHeap || !contribSrcOfs || !numHeaps)
return 0;
pmap -> numContribSrc = 0;
for (heap = 0; heap < numHeaps; heap++) {
contribSrcOfs [heap] = pmap -> numContribSrc;
if (fseek(contribSrcHeap [heap], 0, SEEK_END) < 0)
error(SYSTEM,
"failed photon contribution source seek in buildContribSources()"
);
pmap -> numContribSrc += heapLen = ftell(contribSrcHeap [heap]) /
sizeof(PhotonContribSource);
if (!(pmap -> contribSrc = realloc(pmap -> contribSrc,
pmap -> numContribSrc * sizeof(PhotonContribSource)
)))
error(SYSTEM,
"failed photon contrib source alloc in buildContribSources()"
);
rewind(contribSrcHeap [heap]);
if (fread(pmap -> contribSrc + contribSrcOfs [heap],
sizeof(PhotonContribSource), heapLen, contribSrcHeap [heap]
) != heapLen)
error(SYSTEM,
"failed reading photon contrib source in buildContribSources()"
);
fclose(contribSrcHeap [heap]);
unlink(contribSrcHeapFname [heap]);
}
return pmap -> numContribSrc;
}
static MODCONT *photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Locate photons near ray -> rop which originated from the light source
modifier indexed by ray -> rsrc, and accumulate their contributions
in pmap -> srcContrib. Return total contributions in irrad and pointer
to binned contributions (or NULL if none were found). */
{
const OBJREC *srcMod = findmaterial(source [ray -> rsrc].so);
MODCONT *contrib = (MODCONT*)lu_find(
pmap -> contribTab, srcMod -> oname
) -> data;
Photon *photon;
COLOR flux;
PhotonSearchQueueNode *sqn;
float r2, invArea;
unsigned i;
/* Zero bins for this source modifier */
memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins);
setcolor(irrad, 0, 0, 0);
if (!contrib || !pmap -> maxGather)
/* Modifier not in LUT or zero bandwidth */
return NULL;
/* Lookup photons */
pmap -> squeue.tail = 0;
/* Pass light source index to filter in findPhotons() via
pmap -> lastContribSrc */
pmap -> lastContribSrc.srcIdx = ray -> rsrc;
findPhotons(pmap, ray);
/* Need at least 2 photons */
if (pmap -> squeue.tail < 2) {
#ifdef PMAP_NONEFOUND
sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
ray -> ro ? ray -> ro -> oname : "<null>",
ray -> rop [0], ray -> rop [1], ray -> rop [2]
);
error(WARNING, errmsg);
#endif
return NULL;
}
/* Average radius^2 between furthest two photons to improve accuracy
and get inverse search area 1 / (PI * r^2). */
/* XXX: OMIT extra normalisation factor 1/PI for ambient calculation? */
sqn = pmap -> squeue.node + 1;
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
invArea = 1 / (PI * PI * r2);
/* Skip the extra photon */
for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
/* Get photon's contribution to density estimate */
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
getPhotonFlux(photon, flux);
scalecolor(flux, invArea);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on photon distance */
scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
#endif
addcolor(irrad, flux);
addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux);
}
return contrib;
}
static void preComputeContrib (PhotonMap *pmap)
/* Precompute contributions for a random subset of finalGather *
pmap -> numPhotons photons, and build the precomputed contribution 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;
MODCONT *binnedContribs;
PhotonMap nuPmap;
repComplete = numPreComp = finalGather * pmap -> numPhotons;
if (verbose) {
sprintf(errmsg, "\nPrecomputing contributions for %ld 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++) {
do {
/* 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 and normal at photon
position. Set emitting light source index from origin. */
VCOPY(ray.rop, photon.pos);
ray.rsrc = photonSrcIdx(pmap, &photon);
for (j = 0; j < 3; j++)
ray.ron [j] = photon.norm [j] / 127.0;
/* Get contributions at photon position; retry with another photon
* if no contribs found */
binnedContribs = photonContrib(pmap, &ray, ray.rcol);
} while (!binnedContribs);
/* Append photon to new heap from ray */
newPhoton(&nuPmap, &ray);
/* TODO: Apply wavelet xform to binned contribs, apply lossy compression,
encode to RGBE */
#if 0
encodeContrib(binnedContribs, pmap -> compRatio);
#endif
/* 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 contribution photon map\n");
#if NIX
fflush(stderr);
#endif
}
/* Rebuild underlying data structure, destroying heap */
buildPhotonMap(pmap, NULL, NULL, 1);
/* TODO: Save wavelet coeffs here? */
}
/* Defs for photon emission counter array passed by sub-processes to parent
* via shared memory */
typedef unsigned long PhotonContribCnt;
/* Indices for photon emission counter array: num photons stored and num
* emitted per source */
#define PHOTONCNT_NUMPHOT 0
#define PHOTONCNT_NUMEMIT(n) (1 + n)
/* Photon distribution counter update interval expressed as bitmask;
counters shared among subling subprocesses will only be updated
in multiples of PMAP_CNTUPDATE in order to reduce contention */
#define PMAP_CNTUPDATE 0xffL
void distribPhotonContrib (PhotonMap* pmap, LUTAB *contribTab,
unsigned numContribs, unsigned numProc
)
{
EmissionMap emap;
char errmsg2 [128], shmFname [PMAP_TMPFNLEN];
unsigned srcIdx, proc;
int shmFile, stat, pid;
double *srcFlux, /* Emitted flux per light source */
srcDistribTarget; /* Target photon count
per source */
PhotonContribCnt *photonCnt, /* Photon emission counter array */
lastPhotonCnt [PHOTONCNT_NUMEMIT(nsources)];
unsigned photonCntSize = (sizeof(PhotonContribCnt) *
PHOTONCNT_NUMEMIT(nsources)
);
FILE **contribSrcHeap = NULL;
char **contribSrcHeapFname = NULL;
PhotonContribSourceIdx *contribSrcOfs = NULL;
pid_t procPids [PMAP_MAXPROC];
if (!pmap)
error(USER, "no contribution photon map specified");
if (!nsources)
error(USER, "no light sources");
if (nsources > MAXMODLIST)
error(USER, "too many light sources");
if (!contribTab || !numContribs)
error(USER, "no modifiers specified for contribution photon map");
/* Allocate photon flux per light source; this differs for every
* source as all sources contribute the same number of distributed
* photons (srcDistribTarget), hence the number of photons emitted per
* source does not correlate with its emitted flux. The resulting flux
* per photon is therefore adjusted individually for each source. */
if (!(srcFlux = calloc(nsources, sizeof(double))))
error(SYSTEM, "can't allocate source flux in distribPhotonContrib");
/* ===================================================================
* INITIALISATION - Set up emission and scattering funcs
* =================================================================== */
emap.samples = NULL;
emap.src = NULL;
emap.maxPartitions = MAXSPART;
emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
if (!emap.partitions)
error(USER, "can't allocate source partitions in distribPhotonContrib");
/* Initialise contrib photon map */
initPhotonMap(pmap, PMAP_TYPE_CONTRIB);
initPmapContrib(contribTab);
initPhotonHeap(pmap);
initPhotonEmissionFuncs();
initPhotonScatterFuncs();
/* Per-subprocess / per-source target counts */
pmap -> distribTarget /= numProc;
srcDistribTarget = nsources
? (double)pmap -> distribTarget / nsources
: 0;
if (!pmap -> distribTarget)
error(INTERNAL, "no photons to distribute in distribPhotonContrib");
/* 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, photonCntSize) < 0)
error(SYSTEM, "failed shared mem init in distribPhotonContrib");
photonCnt = mmap(NULL, photonCntSize,
PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0
);
if (photonCnt == MAP_FAILED)
error(SYSTEM, "failed shsared mem mapping in distribPhotonContrib");
#else
/* Allocate photon counters statically on Windoze */
if (!(photonCnt = malloc(photonCntSize)))
error(SYSTEM, "failed trivial malloc in distribPhotonContrib");
for (srcIdx = 0; srcIdx < PHOTONCNT_NUMEMIT(nsources); srcIdx++)
photonCnt [srcIdx] = 0;
#endif /* NIX */
/* Zero overflow tracking counters */
memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt));
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 flux emitted from sources/ports
* ============================================================= */
for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
unsigned portCnt = 0;
const OBJREC *srcMod = findmaterial(source [srcIdx].so);
srcFlux [srcIdx] = 0;
/* Skip this source if its contributions are not sought */
if (!lu_find(pmap -> contribTab, srcMod -> oname) -> data) {
sprintf(errmsg, "ignoring contributions from source %s",
source [srcIdx].so -> oname
);
error(WARNING, errmsg);
continue;
}
emap.src = source + srcIdx;
do { /* Need at least one iteration if no ports! */
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);
srcFlux [srcIdx] += colorAvg(emap.partFlux);
}
portCnt++;
} while (portCnt < numPhotonPorts);
if (srcFlux [srcIdx] < FTINY) {
sprintf(errmsg, "source %s has zero emission",
source [srcIdx].so -> oname
);
error(WARNING, errmsg);
}
}
/* Allocate & init per-subprocess contribution source heap files */
contribSrcHeap = calloc(numProc, sizeof(FILE*));
contribSrcHeapFname = calloc(numProc, sizeof(char*));
contribSrcOfs = calloc(numProc, sizeof(PhotonContribSourceIdx));
if (!contribSrcHeap || !contribSrcHeapFname || !contribSrcOfs)
error(SYSTEM, "failed contribution source heap allocation "
"in distribPhotonContrib()"
);
for (proc = 0; proc < numProc; proc++) {
contribSrcHeapFname [proc] = malloc(PMAP_TMPFNLEN);
if (!contribSrcHeapFname [proc])
error(SYSTEM, "failed contribution source heap file allocation "
"in distribPhotonContrib()"
);
mktemp(strcpy(contribSrcHeapFname [proc], PMAP_TMPFNAME));
if (!(contribSrcHeap [proc] = fopen(contribSrcHeapFname [proc], "w+b")))
error(SYSTEM, "failed opening contribution source heap file "
"in distribPhotonContrib()"
);
}
/* 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; opened and mmapped files inherited */
#else
if (1) {
/* No subprocess under Windoze */
#endif
/* Local photon counters for this subprocess */
unsigned long lastNumPhotons = 0, localNumEmitted = 0;
double photonFluxSum = 0; /* Accum. photon flux */
/* 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 PMAP_SIGUSR
{
double partNumEmit;
unsigned long partEmitCnt;
double srcPhotonFlux, avgPhotonFlux;
unsigned portCnt, passCnt, prePassCnt;
float srcPreDistrib;
double srcNumEmit; /* # to emit from source */
unsigned long srcNumDistrib; /* # stored */
void sigUsrDiags()
/* Loop diags via SIGUSR1 */
{
sprintf(errmsg,
"********************* Proc %d Diags *********************\n"
"srcIdx = %d (%s)\nportCnt = %d (%s)\npassCnt = %d\n"
"srcFlux = %f\nsrcPhotonFlux = %f\navgPhotonFlux = %f\n"
"partNumEmit = %f\npartEmitCnt = %lu\n\n",
proc, srcIdx, findmaterial(source [srcIdx].so) -> oname,
portCnt, photonPorts [portCnt].so -> oname,
passCnt, srcFlux [srcIdx], srcPhotonFlux, avgPhotonFlux,
partNumEmit, partEmitCnt
);
eputs(errmsg);
fflush(stderr);
}
}
#endif
#ifdef PMAP_SIGUSR
signal(SIGUSR1, sigUsrDiags);
#endif
#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
/* =============================================================
* 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
* ============================================================= */
for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
const unsigned numEmitIdx = PHOTONCNT_NUMEMIT(srcIdx);
#ifndef PMAP_SIGUSR
unsigned portCnt, passCnt = 0, prePassCnt = 0;
float srcPreDistrib = preDistrib;
double srcNumEmit = 0; /* # to emit from source */
unsigned long srcNumDistrib = pmap -> numPhotons; /* #stored */
#else
passCnt = prePassCnt = 0;
srcPreDistrib = preDistrib;
srcNumEmit = 0; /* # to emit from source */
srcNumDistrib = pmap -> numPhotons; /* # stored */
#endif
if (srcFlux [srcIdx] < FTINY)
/* Source has zero emission or was skipped in prepass
because its contributions are not sought */
continue;
while (passCnt < 2) {
if (!passCnt) {
/* INIT PASS 1 */
if (++prePassCnt > maxPreDistrib) {
/* Warn if no photons contributed after sufficient
* iterations; only output from subprocess 0 to reduce
* console clutter */
if (!proc) {
sprintf(errmsg,
"source %s: too many prepasses, skipped",
source [srcIdx].so -> oname
);
error(WARNING, errmsg);
}
break;
}
/* Num to emit is fraction of target count */
srcNumEmit = srcPreDistrib * srcDistribTarget;
}
else {
/* INIT PASS 2 */
#ifndef PMAP_SIGUSR
double srcPhotonFlux, avgPhotonFlux;
#endif
/* Based on the outcome of the predistribution we can now
* figure out how many more photons we have to emit from
* the current source to meet the target count,
* srcDistribTarget. This value is clamped to 0 in case
* the target has already been exceeded in pass 1.
* srcNumEmit and srcNumDistrib is the number of photons
* emitted and distributed (stored) from the current
* source in pass 1, respectively. */
srcNumDistrib = pmap -> numPhotons - srcNumDistrib;
srcNumEmit *= srcNumDistrib
? max(srcDistribTarget/srcNumDistrib, 1) - 1
: 0;
if (!srcNumEmit)
/* No photons left to distribute in main pass */
break;
srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit;
avgPhotonFlux = photonFluxSum / (srcIdx + 1);
if (avgPhotonFlux > FTINY &&
srcPhotonFlux / avgPhotonFlux < FTINY
) {
/* Skip source if its photon flux is grossly below the
* running average, indicating negligible contributions
* at the expense of excessive distribution time; only
* output from subproc 0 to reduce console clutter */
if (!proc) {
sprintf(errmsg,
"source %s: itsy bitsy photon flux, skipped",
source [srcIdx].so -> oname
);
error(WARNING, errmsg);
}
srcNumEmit = 0; /* Or just break??? */
}
/* Update sum of photon flux per light source */
photonFluxSum += srcPhotonFlux;
}
portCnt = 0;
do { /* Need at least one iteration if no ports! */
emap.src = source + srcIdx;
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++
) {
#ifndef PMAP_SIGUSR
double partNumEmit;
unsigned long partEmitCnt;
#endif
/* 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;
* scale according to its normalised contribushunn to
* the emitted source flux */
partNumEmit = (srcNumEmit * colorAvg(emap.partFlux) /
srcFlux [srcIdx]
);
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 according to PDF (if any). If
* accepted, allocate associated contribution
* origin, and trace through scene until
* absorbed/leaked; emitPhoton() sets the emitting
* light source index in photonRay */
/* 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;
newPhotonContribSource(pmap, &photonRay,
contribSrcHeap [proc]
);
/* Set subprocess index in photonRay for post-
* distrib contribution source index linearisation;
* this is propagated with the contrib source index
* in photonRay and set for photon hits by
* newPhoton() */
PMAP_SETRAYPROC(&photonRay, proc);
tracePhoton(&photonRay);
/* Update local emission counter */
localNumEmitted++;
if (!(localNumEmitted & PMAP_CNTUPDATE)) {
/* Update global counters shared with siblings
in intervals to reduce overhead/contention */
lastPhotonCnt [numEmitIdx] = photonCnt [numEmitIdx];
photonCnt [numEmitIdx] += PMAP_CNTUPDATE;
/* Check for overflow using previous values */
if (photonCnt [numEmitIdx] <
lastPhotonCnt [numEmitIdx]
) {
sprintf(errmsg, "photon emission counter "
"overflow (was: %ld, is: %ld)",
lastPhotonCnt [numEmitIdx],
photonCnt [numEmitIdx]
);
error(INTERNAL, errmsg);
}
/* Differentially increment photon counter */
lastNumPhotons = pmap -> numPhotons;
photonCnt [PHOTONCNT_NUMPHOT] +=
pmap -> numPhotons - lastNumPhotons;
/* 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 [PHOTONCNT_NUMPHOT] <
lastNumPhotons
) {
sprintf(errmsg, "photon counter overlow "
"(was: %ld, is: %ld)", lastNumPhotons,
photonCnt [PHOTONCNT_NUMPHOT]
);
error(INTERNAL, errmsg);
}
}
}
#if !NIX
/* Synchronous progress report on Windoze */
if (!proc && photonRepTime > 0 &&
time(NULL) >= repLastTime + photonRepTime
) {
unsigned s;
repComplete = pmap -> distribTarget * numProc;
repProgress = photonCnt [PHOTONCNT_NUMPHOT];
for (repEmitted = 0, s = 0; s < nsources; s++)
repEmitted += photonCnt [numEmitIdx];
pmapDistribReport();
}
#endif
}
portCnt++;
} while (portCnt < numPhotonPorts);
if (pmap -> numPhotons == srcNumDistrib) {
/* Double predistrib factor in case no photons were stored
* for this source and redo pass 1 */
srcPreDistrib *= 2;
}
else {
/* Now do pass 2 */
passCnt++;
}
}
}
/* Flush heap buffa one final time to prevent data corruption */
flushPhotonHeap(pmap);
/* Flush last contribution origin to origin heap file */
newPhotonContribSource(pmap, NULL, contribSrcHeap [proc]);
/* Heap files closed automatically on exit
fclose(pmap -> heap);
fclose(orgHeap [proc]); */
#ifdef DEBUG_PMAP
sprintf(
errmsg, "Proc %d total %ld photons\n", proc, pmap -> numPhotons
);
eputs(errmsg);
fflush(stderr);
#endif
#ifdef PMAP_SIGUSR
signal(SIGUSR1, SIG_DFL);
#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 distribPhotonContrib()"
);
else
/* Saves 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 */
repComplete = pmap -> distribTarget * numProc;
repProgress = photonCnt [PHOTONCNT_NUMPHOT];
for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++)
repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
/* Get global photon count from shmem updated by subprocs */
pmap -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT];
if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
pmapDistribReport();
#ifdef SIGCONT
else signal(SIGCONT, pmapDistribReport);
#endif
}
#endif /* NIX */
/* ================================================================
* POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
* ================================================================ */
#ifdef SIGCONT
/* Reset signal handler */
signal(SIGCONT, SIG_DFL);
#endif
free(emap.samples);
if (!pmap -> numPhotons)
error(USER, "empty contribution photon map");
/* Load per-subprocess contribution sources into pmap -> contribSrc */
/* Dumb compilers apparently need the char** cast */
pmap -> numContribSrc = buildContribSources(pmap, contribSrcHeap,
(char**)contribSrcHeapFname, contribSrcOfs, numProc
);
if (!pmap -> numContribSrc)
error(INTERNAL, "no contribution sources in contribution photon map");
/* Set photon flux per source */
/* TODO: HOW/WHERE DO WE HANDLE COEFFICIENT MODE??? */
for (srcIdx = 0; srcIdx < nsources; srcIdx++)
srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
#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("\nBuilding contribution photon map...\n");
#if NIX
fflush(stderr);
#endif
}
/* Build underlying data structure; heap is destroyed */
buildPhotonMap(pmap, srcFlux, contribSrcOfs, numProc);
/* Precompute binned photon contributions */
if (verbose)
eputs("\n");
preComputeContrib(contribPmap);
/* Free per-subprocess origin heap files */
for (proc = 0; proc < numProc; proc++)
free(contribSrcHeapFname [proc]);
free(contribSrcHeapFname);
free(contribSrcHeap);
free(contribSrcOfs);
if (verbose)
eputs("\n");
}
c4science · Help