pmapcontrib.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, May 10, 04:27

pmapcontrib.c
View Options

	#ifndef lint
	static const char RCSid[] = "$Id: pmapcontrib.c,v 2.19 2018/11/08 00:54:07 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.19 2018/11/08 00:54:07 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 *addSourceModifier(
	LUTAB srcContribs, unsigned numSrcContrib, char *mod,
	char binParm, char binVal, int binCnt
	)
	/* Add light source modifier mod to contribution lookup table srcContribs,
	and update numSrcContrib. Return initialised contribution data for this
	modifier. */
	{
	LUENT *lutEntry = lu_find(srcContribs, mod);
	MODCONT *contrib;
	EPNODE *eBinVal;

	if (lutEntry -> data) {
	/* Reject duplicate modifiers */
	sprintf(errmsg, "duplicate light source modifier %s", mod);
	error(USER, errmsg);
	}
	if (*numSrcContrib >= 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 addSourceModifier");
	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;
	++(*numSrcContrib);

	return(contrib);
	}



	void addSourceModfile(
	LUTAB srcContribs, unsigned numSrcContrib, char *modFile,
	char binParm, char binVal, int binCnt
	)
	/* Add light source modifiers from file modFile to contribution lookup table
	* srcContribs, and update numSrcContrib.
	* 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 (*numSrcContrib + 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 */
	addSourceModifier(
	srcContribs, numSrcContrib, mods [i], binParm, binVal, binCnt
	);
	}



	static int originBin (LUTAB srcContribs, const RAY ray)
	/* Map origin 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 \|\| !srcContribs)
	return -1;

	src = &source [ray -> rsrc];
	srcMod = findmaterial(src -> so);
	srcCont = (MODCONT*)lu_find(srcContribs, 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 photonRayBin");
	return -1;
	}

	return bin;
	}



	static PhotonOrigin newPhotonOrigin (
	PhotonMap pmap, const RAY orgRay, FILE *orgHeap
	)
	/* Add origin ray for emitted photon and save light source index, and
	* binned direction for contrib photons. The current origin is stored in
	* pmap -> lastContribOrg. If the previous origin spawned photons (i.e.
	* has srcIdx >= 0), it's appended to orgHeap. If orgRay == NULL, the
	* current origin is still flushed, but no new origin is set.
	* Returns updated origin counter pmap -> numContribOrg. */
	{
	if (!pmap \|\| !orgHeap)
	return 0;

	/* Check if last origin ray has spawned photons (srcIdx >= 0, see
	* newPhoton()), in which case we save it to the origin heap file
	* before clobbering it. (Note this is short-term storage, so we doan'
	* need da portable I/O stuff here). */
	if (pmap -> lastContribOrg.srcIdx >= 0) {
	if (!fwrite(
	&pmap -> lastContribOrg, sizeof(PhotonContribOrigin), 1, orgHeap
	))
	error(
	SYSTEM,
	"failed writing contrib photon origin in newPhotonOrigin"
	);

	pmap -> numContribOrg++;
	if (!pmap -> numContribOrg \|\| pmap -> numContribOrg > PMAP_MAXORIGIN)
	/* numContribOrg overflowed (wrapped around) or exceeded max */
	error(INTERNAL, "contrib photon origin overflow in newPhotonOrigin");
	}

	/* Mark this origin unused with a negative source index until its path
	spawns a photon (see newPhoton()) */
	pmap -> lastContribOrg.srcIdx = -1;
	/* Map ray to bin in anticipation that this origin will be used, since the
	ray will be lost once a photon is spawned */
	pmap -> lastContribOrg.srcBin = originBin(pmap -> srcContrib, orgRay);

	return pmap -> numContribOrg;
	}



	static PhotonOrigin buildOrigins (
	PhotonMap pmap, FILE orgHeap, char *orgHeapFname,
	PhotonOrigin *orgOfs, unsigned numHeaps
	)
	/* Consolidate per-subprocess photon origin heaps into the origin array
	* pmap -> contribOrg. Returns offset for origin index linearisation in
	* pmap -> numContribOrg. The heap files in orgHeap are closed on return. */
	{
	PhotonOrigin heapLen;
	unsigned heap;

	if (!pmap \|\| !orgHeap \|\| !orgOfs \|\| !numHeaps)
	return 0;

	pmap -> numContribOrg = 0;

	for (heap = 0; heap < numHeaps; heap++) {
	orgOfs [heap] = pmap -> numContribOrg;

	if (fseek(orgHeap [heap], 0, SEEK_END) < 0)
	error(SYSTEM, "failed contrib photon origin seek in buildOrigins");
	pmap -> numContribOrg +=
	heapLen = ftell(orgHeap [heap]) / sizeof(PhotonContribOrigin);

	if (!(pmap -> contribOrg = realloc(
	pmap -> contribOrg,
	pmap -> numContribOrg * sizeof(PhotonContribOrigin)
	)))
	error(SYSTEM, "failed contrib photon origin alloc in buildOrigins");

	rewind(orgHeap [heap]);
	if (fread(
	pmap -> contribOrg + orgOfs [heap], sizeof(PhotonContribOrigin),
	heapLen, orgHeap [heap]
	) != heapLen)
	error(SYSTEM, "failed reading contrib photon origin in buildOrigins");

	fclose(orgHeap [heap]);
	unlink(orgHeapFname [heap]);
	}

	return pmap -> numContribOrg;
	}



	void 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. */
	{
	const OBJREC *srcMod = findmaterial(source [ray -> rsrc].so);
	MODCONT srcContrib = (MODCONT)lu_find(
	pmap -> srcContrib, srcMod -> oname
	) -> data;
	Photon *photon;
	COLOR flux;
	PhotonSearchQueueNode *sqn;
	float r2, invArea;
	unsigned i;

	/* Zero bins for this source modifier */
	memset(srcContrib -> cbin, 0, sizeof(DCOLOR) * srcContrib -> nbins);
	setcolor(irrad, 0, 0, 0);

	if (!srcContrib \|\| !pmap -> maxGather)
	/* Modifier not in LUT or zero bandwidth */
	return;

	/* Lookup photons */
	pmap -> squeue.tail = 0;
	/* Pass light source index to filter in findPhotons() via
	pmap -> lastContribOrg */
	pmap -> lastContribOrg.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;
	}

	/* 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(srcContrib -> cbin [photonSrcBin(pmap, photon)], flux);
	}

	return;
	}



	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;
	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++) {
	/* 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 */
	photonContrib(pmap, &ray, ray.rcol);

	/* Append photon to new heap from ray */
	newPhoton(&nuPmap, &ray);

	/* TODO: wavelet xform on contribs here? */
	/* waveletEncode(...); */


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



	void distribPhotonContrib (
	PhotonMap* pmap, LUTAB *srcContrib, unsigned numSrcContrib,
	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 */
	unsigned photonCntSize = sizeof(PhotonContribCnt) *
	PHOTONCNT_NUMEMIT(nsources);
	FILE **orgHeap = NULL;
	char **orgHeapFname = NULL;
	PhotonOrigin *orgOfs = NULL;

	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 (!srcContrib \|\| !numSrcContrib)
	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(srcContrib);
	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 shared 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 */

	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 -> srcContrib, 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 origin heap files */
	orgHeap = calloc(numProc, sizeof(FILE*));
	orgHeapFname = calloc(numProc, sizeof(char*));
	orgOfs = calloc(numProc, sizeof(PhotonOrigin));
	if (!orgHeap \|\| !orgHeapFname \|\| !orgOfs)
	error(
	SYSTEM, "failed contribution origin heap allocation "
	"in distribPhotonContrib"
	);

	for (proc = 0; proc < numProc; proc++) {
	orgHeapFname [proc] = malloc(PMAP_TMPFNLEN);
	if (!orgHeapFname [proc])
	error(
	SYSTEM, "failed contribution origin heap file allocation "
	"in distribPhotonContrib"
	);

	mktemp(strcpy(orgHeapFname [proc], PMAP_TMPFNAME));
	if (!(orgHeap [proc] = fopen(orgHeapFname [proc], "w+b")))
	error(
	SYSTEM, "failed opening contribution origin 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++) {
	#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++;

	/* Update local and shared global emission counter */
	photonCnt [PHOTONCNT_NUMEMIT(srcIdx)] += partEmitCnt;
	localNumEmitted += partEmitCnt;

	/* Integer counter avoids FP rounding errors during
	* iteration */
	while (partEmitCnt--) {
	RAY photonRay;

	/* Emit photon according to PDF (if any), allocate
	* associated contribution origin, and trace through
	* scene until absorbed/leaked; emitPhoton() sets the
	* emitting light source index in photonRay */
	emitPhoton(&emap, &photonRay);
	#if 1
	if (emap.port)
	/* !!! PHOTON PORT REJECTION SAMPLING HACK: set
	* !!! photon port as fake hit object for
	* !!! photon ray to check for intersection in
	* !!! tracePhoton() */
	photonRay.ro = emap.port -> so;
	#endif
	newPhotonOrigin(pmap, &photonRay, orgHeap [proc]);
	/* Set subprocess index in photonRay for post-
	* distrib contribution origin index linearisation;
	* this is propagated with the origin index in
	* photonRay and set for photon hits by newPhoton() */
	PMAP_SETRAYPROC(&photonRay, proc);
	tracePhoton(&photonRay);
	}

	/* Update shared global photon count */
	photonCnt [PHOTONCNT_NUMPHOT] +=
	pmap -> numPhotons - lastNumPhotons;
	lastNumPhotons = pmap -> numPhotons;
	#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 [PHOTONCNT_NUMEMIT(s)];

	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 */
	newPhotonOrigin(pmap, NULL, orgHeap [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 if (pid < 0)
	error(SYSTEM, "failed to fork subprocess in distribPhotonContrib");
	}

	#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))
	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 origins into pmap -> contribOrg */
	/* Dumb compilers apparently need the char** cast */
	pmap -> numContribOrg = buildOrigins(
	pmap, orgHeap, (char**)orgHeapFname, orgOfs, numProc
	);
	if (!pmap -> numContribOrg)
	error(INTERNAL, "no origins 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, orgOfs, numProc);

	/* Precompute binned photon contributions */
	if (verbose)
	eputs("\n");
	preComputeContrib(contribPmap);

	/* Free per-subprocess origin heap files */
	for (proc = 0; proc < numProc; proc++)
	free(orgHeapFname [proc]);

	free(orgHeapFname);
	free(orgHeap);
	free(orgOfs);

	if (verbose)
	eputs("\n");
	}

pmapcontrib.cNo OneTemporaryActions

File Metadata

pmapcontrib.cView Options

Event Timeline

pmapcontrib.c
No OneTemporary
Actions

pmapcontrib.c
View Options