pmapcontrib.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, May 12, 05:45

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 "pmapio.h"
	#include "pmapdiag.h"
	#include "otypes.h"
	#include "otspecial.h"
	#if NIX
	#include <sys/mman.h>
	#include <sys/wait.h>
	#endif



	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 int contribSourceBin (LUTAB contribs, const RAY ray)
	/* Map contribution source ray to its bin for light source ray -> rsrc.
	Return invalid bin -1 if mapping failed. */
	{
	const SRCREC *src;
	const OBJREC *srcMod;
	const MODCONT *srcCont;
	double binReal;
	int bin;

	/* 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;

	if (srcCont -> binv -> type != NUM) {
	/* Binning function defined; set up shadow ray pointing to source
	for evaluation */
	RAY srcRay;
	int i;

	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);
	}

	if ((binReal = evalue(srcCont -> binv)) < -.5)
	return -1;

	if ((bin = binReal + .5) >= srcCont -> nbins) {
	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);

	#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");
	}

pmapcontrib.cNo OneTemporaryActions

File Metadata

pmapcontrib.cView Options

Event Timeline

pmapcontrib.c
No OneTemporary
Actions

pmapcontrib.c
View Options