pmap.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Tue, Apr 16, 10:27

pmap.c
View Options

	#ifndef lint
	static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $";
	#endif


	/*
	======================================================================
	Photon map main module

	Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
	(c) Fraunhofer Institute for Solar Energy Systems,
	supported by the German Research Foundation
	(DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS")
	(c) Lucerne University of Applied Sciences and Arts,
	supported by the Swiss National Science Foundation
	(SNSF #147053, "Daylight Redirecting Components")
	(c) Tokyo University of Science,
	supported by the JSPS Grants-in-Aid for Scientific Research
	(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
	======================================================================

	$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $
	*/


	#include "pmap.h"
	#include "pmapmat.h"
	#include "pmapsrc.h"
	#include "pmaprand.h"
	#include "pmaproi.h"
	#include "pmapio.h"
	#include "pmapdens.h"
	#include "pmapbias.h"
	#include "pmapdiag.h"
	#include "pmutil.h"
	#include "otypes.h"
	#include "otspecial.h"
	#include <time.h>
	#if NIX
	#include <sys/stat.h>
	#include <sys/mman.h>
	#include <sys/wait.h>
	#include <sys/file.h>
	#endif



	static int photonParticipate (RAY *ray)
	/* Trace photon through participating medium (including dielectric
	and interface). Returns 1 if passed through, or 0 if absorbed.
	Analogon to rayparticipate(). */
	{
	int i;
	RREAL xi1, cosTheta, phi, du, dv, rmax0;
	const RREAL cext = colorAvg(ray -> cext), invCext = 1 / cext,
	albedo = colorAvg(ray -> albedo), invAlbedo = 1 / albedo,
	gecc = ray -> gecc, gecc2 = sqr(gecc);
	FVECT u, v;
	COLOR cvext;

	/* Save incoming ray's max distance; this is required since rmax is
	set to the mean free distance for each path segment (see below) */
	rmax0 = ray -> rmax;

	/* Mean free distance until interaction with medium */
	ray -> rmax = -log(pmapRandom(mediumState)) / cext;

	while (!localhit(ray, &thescene)) {
	/* No geometry hit within mean free dist ray->rmax; photon interacts
	with medium */
	if (!incube(&thescene, ray -> rop)) {
	/* Terminate photon if it has leaked from the scene */
	#ifdef DEBUG_PMAP_LEAKED
	fprintf(
	stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n",
	ray -> rop [0], ray -> rop [1], ray -> rop [2]
	);
	#endif
	return 0;
	}

	/* RGB attenuation over mean free distance */
	setcolor(cvext,
	exp(-ray -> rmax * ray -> cext [0]),
	exp(-ray -> rmax * ray -> cext [1]),
	exp(-ray -> rmax * ray -> cext [2])
	);

	/* Modify ray color and normalise */
	multcolor(ray -> rcol, cvext);
	colorNorm(ray -> rcol);
	/* Set ray origin for next interaction */
	VCOPY(ray -> rorg, ray -> rop);

	/* Update ray's path length; terminate photon if maximum reached.
	See also photonRay(). */
	rmax0 -= ray -> rot;
	if (photonMaxPathDist > 0 && rmax0 < 0)
	return 0;

	/* Store volume photon if nonzero albedo (or unconditionally if
	lightflow photon); this also accounts for direct inscattering
	from light sources */
	if (albedo > FTINY \|\| lightFlowPhotonMapping) {
	/* Set ray's path length to store with photon */
	ray -> rmax = rmax0;

	if (lightFlowPhotonMapping) {
	/* Temporarily correct normalised flux for lightflow photons
	based on extinction (= mean num photons deposited per unit
	path length). */
	scalecolor(ray -> rcol, invCext);
	newPhoton(lightFlowPmap, ray, NULL);
	/* Undo correction for next iteration after storing photon */
	scalecolor(ray -> rcol, cext);
	}
	else
	newPhoton(volumePmap, ray, NULL);
	}

	/* Colour bleeding without attenuation (?) */
	multcolor(ray -> rcol, ray -> albedo);
	scalecolor(ray -> rcol, invAlbedo);

	if (pmapRandom(rouletteState) < albedo && !lightFlowPhotonMapping) {
	/* Scatter volume photon; lightflow photons are NOT scattered!

	NOTE: "albedo" here specifically refers to the SCATTERING ALBEDO
	sigma_s / sigma_t (see Greg's email "Albedo within dielectric /
	interface material?" dated 24 Jun 2021), where sigma_s is the
	is the scattering probability, and sigma_t is the extinction.
	The extinction sigma_t = sigma_a + sigma_s is itself the
	sum of the absorption and scattering probabilities, sigma_a
	resp. sigma_s. */
	xi1 = pmapRandom(scatterState);
	cosTheta = ray -> gecc <= FTINY
	? 2 * xi1 - 1
	: 0.5 / gecc * (
	1 + gecc2 - sqr((1 - gecc2) / (1 + gecc * (2 * xi1 - 1)))
	);

	phi = 2 * PI * pmapRandom(scatterState);
	du = dv = sqrt(1 - sqr(cosTheta)); /* sin(theta) */
	du *= cos(phi);
	dv *= sin(phi);

	/* Get axes u & v perpendicular to photon direction */
	i = 0;
	do {
	v [0] = v [1] = v [2] = 0;
	v [i++] = 1;
	fcross(u, v, ray -> rdir);
	} while (normalize(u) < FTINY);
	fcross(v, ray -> rdir, u);

	for (i = 0; i < 3; i++)
	ray -> rdir [i] = du * u [i] + dv * v [i] +
	cosTheta * ray -> rdir [i];
	}
	#if 0
	/* Photons should be terminated if not scattered; however, this
	unconditionally terminates all photon in dielectrics, since their
	scattering albedo is implicitly zero */
	else return 0;
	#endif

	ray -> rlvl++;
	ray -> rmax = -log(pmapRandom(mediumState)) / cext;
	/* fprintf(stderr, "%.3f\n", ray -> rmax); */
	}

	/* Passed through medium until intersecting local object */
	setcolor(cvext,
	exp(-ray -> rot * ray -> cext [0]),
	exp(-ray -> rot * ray -> cext [1]),
	exp(-ray -> rot * ray -> cext [2])
	);

	/* Modify ray color and normalise */
	multcolor(ray -> rcol, cvext);
	colorNorm(ray -> rcol);

	/* Restore outgoing ray's max distance */
	ray -> rmax = rmax0;

	return 1;
	}



	void tracePhoton (RAY *ray)
	/* Follow photon as it bounces around... */
	{
	long mod;
	OBJREC mat, port = NULL;

	if (!ray -> parent) {
	/* !!! PHOTON PORT REJECTION SAMPLING HACK: get photon port for
	* !!! primary ray from ray -> ro, then reset the latter to NULL so
	* !!! as not to interfere with localhit() */
	port = ray -> ro;
	ray -> ro = NULL;
	}

	if (ray -> rlvl > photonMaxBounce) {
	#ifdef PMAP_RUNAWAY_WARN
	error(WARNING, "runaway photon!");
	#endif
	return;
	}

	if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
	return;

	/* NOTE: localhit() limits intersections to the aft plane's distance
	ray->rmax. This mechanism is (ab)used here to limit the photon path
	length by initialising ray->rmax to photonMaxPathDist; see emitPhoton().
	With unlimited path length (photonMaxPathDist=0), ray->rmax becomes
	negative (see photonRay()), and the aft plane has no effect. */
	if (localhit(ray, &thescene)) {
	mod = ray -> ro -> omod;
	if (port && ray -> ro != port) {
	/* Terminate photon if emitted from port without intersecting it */
	#ifdef PMAP_PORTREJECT_WARN
	sprintf(errmsg, "photon outside port %s", ray -> ro -> oname);
	error(WARNING, errmsg);
	#endif
	return;
	}

	if ((ray -> clipset && inset(ray -> clipset, mod)) \|\| mod == OVOID) {
	/* Transfer ray if modifier is VOID or clipped within antimatta */
	RAY tray;
	photonRay(ray, &tray, PMAP_XFER, NULL);
	tracePhoton(&tray);
	}
	else {
	/* Scatter for modifier material */
	mat = objptr(mod);
	photonScatter [mat -> otype] (mat, ray);
	}
	}
	}



	static void preComputeGlobal (PhotonMap *pmap)
	/* Precompute irradiance from global photons for final gathering for
	a random subset of finalGather * pmap -> numPhotons photons, and builds
	the photon map, discarding the original photons. */
	/* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */
	{
	unsigned long p, numPreComp;
	unsigned j;
	const double pInc = finalGather > FTINY ? 1 / finalGather : 0;
	PhotonIdx pIdx;
	Photon photon;
	RAY ray;
	PhotonMap nuPmap;

	repComplete = numPreComp = finalGather * pmap -> numPhotons;

	if (verbose) {
	sprintf(errmsg,
	"\nPrecomputing irradiance for %ld global photons\n", numPreComp
	);
	eputs(errmsg);
	#if NIX
	fflush(stderr);
	#endif
	}

	/* Copy photon map for precomputed photons */
	memcpy(&nuPmap, pmap, sizeof(PhotonMap));

	/* Zero counters, init new heap and extents */
	nuPmap.numPhotons = 0;
	initPhotonHeap(&nuPmap);

	for (j = 0; j < 3; j++) {
	nuPmap.minPos [j] = FHUGE;
	nuPmap.maxPos [j] = -FHUGE;
	}

	/* Record start time, baby */
	repStartTime = repLastTime = time(NULL);
	#ifdef SIGCONT
	signal(SIGCONT, pmapPreCompReport);
	#endif
	repProgress = 0;

	photonRay(NULL, &ray, PRIMARY, NULL);
	ray.ro = NULL;

	for (p = 0; p < numPreComp; p++) {
	/* Get random photon from stratified distribution in source heap to
	* avoid duplicates and clustering.
	* NOTE: multiplication with pInc expanded to prevent rounding error
	* which would occasionally cause the last photon to be selected
	* for precomputation again! */
	pIdx = firstPhoton(pmap) + (unsigned long)(
	p * pInc + pmapRandom(pmap -> randState) * pInc
	);
	getPhoton(pmap, pIdx, &photon);

	/* Init dummy photon ray with intersection at photon position */
	VCOPY(ray.rop, photon.pos);
	for (j = 0; j < 3; j++)
	ray.ron [j] = photon.norm [j] / 127.0;

	/* Get density estimate at photon position */
	photonDensity(pmap, &ray, ray.rcol);

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

	/* Update progress */
	repProgress++;

	if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
	pmapPreCompReport();
	#ifdef SIGCONT
	else signal(SIGCONT, pmapPreCompReport);
	#endif
	}

	/* Flush heap */
	flushPhotonHeap(&nuPmap);

	#ifdef SIGCONT
	signal(SIGCONT, SIG_DFL);
	#endif

	/* Trash original pmap, replace with precomputed one */
	deletePhotons(pmap);
	memcpy(pmap, &nuPmap, sizeof(PhotonMap));

	if (verbose) {
	eputs("\nRebuilding precomputed photon map\n");
	#if NIX
	fflush(stderr);
	#endif
	}

	/* Rebuild underlying data structure, destroying heap */
	buildPhotonMap(pmap, NULL, NULL, 1);
	}



	typedef struct {
	unsigned long numPhotons [NUM_PMAP_TYPES],
	numEmitted, numComplete;
	} PhotonCnt;



	void distribPhotons (PhotonMap **pmaps, unsigned numProc)
	{
	EmissionMap emap;
	char errmsg2 [128], shmFname [PMAP_TMPFNLEN];
	unsigned t, srcIdx, proc;
	double totalFlux = 0;
	int shmFile, stat, pid;
	PhotonMap *pm;
	PhotonCnt *photonCnt;
	pid_t procPids [PMAP_MAXPROC];

	for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++);

	if (t >= NUM_PMAP_TYPES)
	error(USER, "no photon maps defined in distribPhotons");

	if (!nsources)
	error(USER, "no light sources in distribPhotons");

	/* ===================================================================
	* INITIALISATION - Set up emission and scattering funcs
	* =================================================================== */
	emap.samples = NULL;
	emap.maxPartitions = MAXSPART;
	emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
	if (!emap.partitions)
	error(INTERNAL, "can't allocate source partitions in distribPhotons");

	/* Initialise all defined photon maps */
	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if (pmaps [t]) {
	initPhotonMap(pmaps [t], t);
	/* Open photon heapfile */
	initPhotonHeap(pmaps [t]);
	/* Per-subprocess target count */
	pmaps [t] -> distribTarget /= numProc;

	if (!pmaps [t] -> distribTarget)
	error(INTERNAL, "no photons to distribute in distribPhotons");
	}

	initPhotonEmissionFuncs();
	initPhotonScatterFuncs();

	/* Get photon ports from modifier list */
	getPhotonPorts(photonPortList);

	/* Get photon sensor modifiers */
	getPhotonSensors(photonSensorList);

	/* Get polyhedral regions of interest */
	getPolyROIs(pmapROImodList);

	#if NIX
	/* Set up shared mem for photon counters (zeroed by ftruncate) */
	strcpy(shmFname, PMAP_TMPFNAME);
	shmFile = mkstemp(shmFname);

	if (shmFile < 0 \|\| ftruncate(shmFile, sizeof(*photonCnt)) < 0)
	error(SYSTEM, "failed shared mem init in distribPhotons");

	photonCnt = mmap(NULL, sizeof(*photonCnt),
	PROT_READ \| PROT_WRITE, MAP_SHARED, shmFile, 0
	);

	if (photonCnt == MAP_FAILED)
	error(SYSTEM, "failed mapping shared memory in distribPhotons");
	#else
	/* Allocate photon counters statically on Windoze */
	if (!(photonCnt = malloc(sizeof(PhotonCnt))))
	error(SYSTEM, "failed trivial malloc in distribPhotons");
	photonCnt -> numEmitted = photonCnt -> numComplete = 0;
	#endif /* NIX */

	if (verbose) {
	sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);

	if (photonPorts) {
	sprintf(errmsg2, " via %d ports", numPhotonPorts);
	strcat(errmsg, errmsg2);
	}

	strcat(errmsg, "\n");
	eputs(errmsg);
	}

	/* ===================================================================
	* FLUX INTEGRATION - Get total photon flux from light sources
	* =================================================================== */
	for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
	unsigned portCnt = 0;

	emap.src = source + srcIdx;
	do {
	/* Set next photon port if defined; note this loop iterates once if
	* no ports are defined. */
	emap.port = emap.src -> sflags & SDISTANT
	? photonPorts + portCnt
	: NULL;
	photonPartition [emap.src -> so -> otype] (&emap);

	if (verbose) {
	sprintf(errmsg, "\tIntegrating flux from source %s ",
	source [srcIdx].so -> oname
	);

	if (emap.port) {
	sprintf(errmsg2, "via port %s ",
	photonPorts [portCnt].so -> oname
	);
	strcat(errmsg, errmsg2);
	}

	sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
	strcat(errmsg, errmsg2);
	eputs(errmsg);
	#if NIX
	fflush(stderr);
	#endif
	}

	for (emap.partitionCnt = 0;
	emap.partitionCnt < emap.numPartitions;
	emap.partitionCnt++
	) {
	initPhotonEmission(&emap, pdfSamples);
	totalFlux += colorAvg(emap.partFlux);
	}

	portCnt++;
	} while (portCnt < numPhotonPorts);
	}

	if (totalFlux < FTINY)
	error(USER, "zero flux from light sources");

	/* Record start time for progress reports */
	repStartTime = time(NULL);

	if (verbose) {
	sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
	eputs(errmsg);
	}

	/* MAIN LOOP */
	for (proc = 0; proc < numProc; proc++) {
	#if NIX
	if (!(pid = fork())) {
	/* SUBPROCESS ENTERS HERE; open and mmapped files inherited */
	#else
	if (1) {
	/* No subprocess under Windoze */
	#endif
	/* Local photon counters for this subprocess */
	unsigned passCnt = 0, prePassCnt = 0;
	unsigned long preIncPhotonCnt, lastNumPhotons [NUM_PMAP_TYPES],
	numEmitted = 0, lastNumEmitted = 0;
	/* Decorrelate shared photon counter updates to reduce
	* contention by setting unique update intervals per subproc */
	const unsigned photonCntUpdate = PMAP_CNTUPDATE + proc;

	/* Seed RNGs from PID for decorellated photon distribution */
	pmapSeed(randSeed + proc, partState);
	pmapSeed(randSeed + (proc + 1) % numProc, emitState);
	pmapSeed(randSeed + (proc + 2) % numProc, cntState);
	pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
	pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
	pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);

	#ifdef DEBUG_PMAP
	/* Output child process PID after random delay to prevent corrupted
	* console output due to race condition */
	usleep(1e6 * pmapRandom(rouletteState));
	fprintf(stderr,
	"Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n",
	proc, getpid()
	);
	/* Allow time for debugger to attach to child process */
	sleep(10);
	#endif

	for (t = 0; t < NUM_PMAP_TYPES; t++)
	lastNumPhotons [t] = 0;

	/* =============================================================
	* 2-PASS PHOTON DISTRIBUTION
	* Pass 1 (pre): emit fraction of target photon count
	* Pass 2 (main): based on outcome of pass 1, estimate remaining
	* number of photons to emit to approximate target
	* count
	* ============================================================= */
	do {
	double numEmit;

	if (!passCnt) {
	/* INIT PASS 1 */
	/* Skip if no photons contributed after sufficient
	* iterations; make it clear to user which photon maps are
	* missing so (s)he can check geometry and materials */
	if (++prePassCnt > maxPreDistrib) {
	sprintf(errmsg, "proc %d: too many prepasses", proc);

	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if (pmaps [t] && !pmaps [t] -> numPhotons) {
	sprintf(errmsg2, ", no %s photons stored",
	pmapName [t]
	);
	strcat(errmsg, errmsg2);
	}

	sprintf(errmsg2, "; increase -apD");
	strcat(errmsg, errmsg2);
	error(USER, errmsg);
	break;
	}

	/* Num to emit is fraction of minimum target count */
	numEmit = FHUGE;

	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if (pmaps [t])
	numEmit = min(pmaps [t] -> distribTarget, numEmit);

	numEmit *= preDistrib;
	}
	else {
	/* INIT PASS 2 */
	/* Based on the outcome of the predistribution we can now
	* estimate how many more photons we have to emit for each
	* photon map to meet its respective target count. This
	* value is clamped to 0 in case the target has already been
	* exceeded in the pass 1. */
	double maxDistribRatio = 0;

	/* Set the distribution ratio for each map; this indicates
	* how many photons of each respective type are stored per
	* emitted photon, and is used as probability for storing a
	* photon by newPhoton(). Since this biases the photon
	* density, newPhoton() promotes the flux of stored photons
	* to compensate. */
	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if ((pm = pmaps [t])) {
	pm -> distribRatio =
	(double)pm -> distribTarget / pm -> numPhotons - 1;

	/* Check if photon map "overflowed", i.e. exceeded its
	* target count in the prepass; correcting the photon
	* flux via the distribution ratio is no longer
	* possible, as no more photons of this type will be
	* stored, so notify the user rather than deliver
	* incorrect results. In future we should handle this
	* more intelligently by using the photonFlux in each
	* photon map to individually correct the flux after
	* distribution. */
	if (pm -> distribRatio <= FTINY) {
	sprintf(errmsg,
	"%s photon map overflow in prepass, reduce -apD ",
	pmapName [t]
	);
	error(INTERNAL, errmsg);
	}

	maxDistribRatio = max(pm -> distribRatio,
	maxDistribRatio
	);
	}

	/* Normalise distribution ratios and calculate number of
	* photons to emit in main pass */
	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if ((pm = pmaps [t]))
	pm -> distribRatio /= maxDistribRatio;

	if ((numEmit = numEmitted * maxDistribRatio) < FTINY)
	/* No photons left to distribute in main pass */
	break;
	}

	/* PHOTON DISTRIBUTION LOOP */
	for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
	unsigned portCnt = 0;
	emap.src = source + srcIdx;

	do {
	/* Set next photon port if defined; note this loop iterates
	* once if no ports are defined. */
	emap.port = emap.src -> sflags & SDISTANT
	? photonPorts + portCnt
	: NULL;
	photonPartition [emap.src -> so -> otype] (&emap);

	if (verbose && !proc) {
	/* Output from subproc 0 only to avoid race condition
	* on console I/O */
	if (!passCnt)
	sprintf(errmsg, "\tPREPASS %d on source %s ",
	prePassCnt, source [srcIdx].so -> oname
	);
	else
	sprintf(errmsg, "\tMAIN PASS on source %s ",
	source [srcIdx].so -> oname
	);

	if (emap.port) {
	sprintf(errmsg2, "via port %s ",
	photonPorts [portCnt].so -> oname
	);
	strcat(errmsg, errmsg2);
	}

	sprintf(errmsg2, "(%lu partitions)\n",
	emap.numPartitions
	);
	strcat(errmsg, errmsg2);
	eputs(errmsg);
	#if NIX
	fflush(stderr);
	#endif
	}

	for (emap.partitionCnt = 0;
	emap.partitionCnt < emap.numPartitions;
	emap.partitionCnt++
	) {
	double partNumEmit;
	unsigned long partEmitCnt;

	/* Get photon origin within current source partishunn
	* and build emission map */
	photonOrigin [emap.src -> so -> otype] (&emap);
	initPhotonEmission(&emap, pdfSamples);

	/* Number of photons to emit from ziss partishunn --
	* proportional to flux; photon ray weight and scalar
	* flux are uniform (latter only varying in RGB). */
	partNumEmit = numEmit * colorAvg(emap.partFlux) /
	totalFlux;
	partEmitCnt = (unsigned long)partNumEmit;

	/* Probabilistically account for fractional photons */
	if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
	partEmitCnt++;

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

	/* Emit photon based on PDF, skip if rejected. */
	/* NOTE: rejection sampling skips incrementing the
	* emission counter (see below), thus compensating
	* for the rejected photons by increasing the photon
	* flux in proportion to the lower effective
	* emission count.
	* BUG: THIS INTERFERES WITH THE PROGRESS COUNTER
	* REPORTED TO THE PARENT, AND WITH THE
	* PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL
	* FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED!
	*/
	if (!emitPhoton(&emap, &photonRay))
	continue;

	/* Set subprocess ID in photonRay; extends path ID.
	Used by newPhoton() as photon attributes. */
	PMAP_SETRAYPROC(&photonRay, proc);
	tracePhoton(&photonRay);

	/* Update local emission counter*/
	numEmitted++;

	if (!(partEmitCnt && (numEmitted & photonCntUpdate))) {
	/* Update global counters shared with siblings in
	intervals to reduce overhead & contention;
	ALSO WHEN EMITTING FINAL PHOTON FROM THIS
	PARTITION, OTHERWISE COUNTERS WILL SKIP LAST
	UPDATE, BIASING THE PHOTON FLUX! */
	#if NIX
	/* Shared photon counter file must be locked! */
	shmLock(shmFile, F_WRLCK);
	#endif
	/* Differentially increment number emitted using
	* local counter numEmitted */
	preIncPhotonCnt = photonCnt -> numEmitted;
	photonCnt -> numEmitted += numEmitted - lastNumEmitted;
	photonCnt -> numComplete += numEmitted - lastNumEmitted;
	lastNumEmitted = numEmitted;

	/* Check overflow using pre-increment values */
	if (photonCnt -> numEmitted < preIncPhotonCnt) {
	sprintf(errmsg, "photon emission counter "
	"overflow (was: %ld, is: %ld)",
	preIncPhotonCnt, photonCnt -> numEmitted
	);
	error(INTERNAL, errmsg);
	}

	/* Update global photon count for each pmap */
	for (t = 0; t < NUM_PMAP_TYPES; t++) {
	if (pmaps [t]) {
	preIncPhotonCnt = photonCnt->numPhotons [t];
	/* Differentially increment using local
	counter lastNumPhotons */
	photonCnt -> numPhotons [t] +=
	pmaps [t] -> numPhotons - lastNumPhotons [t];
	lastNumPhotons [t] = pmaps [t] -> numPhotons;

	/* Check for photon counter overflow (this
	* could only happen before an emission
	* counter overflow if the scene has an
	* absurdly high albedo and/or very dense
	* geometry) */
	if (photonCnt -> numPhotons [t] < preIncPhotonCnt) {
	sprintf(errmsg, "photon counter "
	"overflow (was: %ld, is: %ld)",
	preIncPhotonCnt, photonCnt -> numPhotons [t]
	);
	error(INTERNAL, errmsg);
	}
	}
	}
	#if NIX
	/* Release lock on shared photon counter file! */
	shmLock(shmFile, F_UNLCK);
	#endif
	}
	}
	#if !NIX
	/* Synchronous progress report on Windoze */
	if (!proc && photonRepTime > 0 &&
	time(NULL) >= repLastTime + photonRepTime) {
	repEmitted = repProgress = photonCnt -> numEmitted;
	repComplete = photonCnt -> numComplete;

	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if (pmaps [t]) {
	repProgress += pmaps [t] -> numPhotons;
	repComplete += pmaps [t] -> distribTarget;
	}
	pmapDistribReport();
	}
	#endif
	}

	portCnt++;
	} while (portCnt < numPhotonPorts);
	}

	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if (pmaps [t] && !pmaps [t] -> numPhotons) {
	/* Double preDistrib in case a photon map is empty and
	* redo pass 1 --> possibility of infinite loop for
	* pathological scenes (e.g. absorbing materials) */
	preDistrib *= 2;
	break;
	}

	if (t >= NUM_PMAP_TYPES)
	/* No empty photon maps found; now do pass 2 */
	passCnt++;
	} while (passCnt < 2);

	/* Flush heap buffa for every pmap one final time;
	* avoids potential data corruption! */
	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if (pmaps [t]) {
	flushPhotonHeap(pmaps [t]);
	/* Heap file closed automatically on exit
	fclose(pmaps [t] -> heap); */
	#ifdef DEBUG_PMAP
	sprintf(errmsg, "Proc %d: total %ld photons\n",
	proc, pmaps [t] -> numPhotons
	);
	eputs(errmsg);
	#endif
	}
	#if NIX
	/* Terminate subprocess */
	exit(0);
	#endif
	}
	else {
	/* PARENT PROCESS CONTINUES LOOP ITERATION HERE */
	if (pid < 0)
	error(SYSTEM, "failed to fork subprocess in distribPhotons");
	else
	/* Save child process IDs */
	procPids [proc] = pid;
	}
	}

	#if NIX
	/* PARENT PROCESS CONTINUES HERE */
	#ifdef SIGCONT
	/* Enable progress report signal handler */
	signal(SIGCONT, pmapDistribReport);
	#endif
	/* Wait for subprocesses to complete while reporting progress */
	proc = numProc;
	while (proc) {
	while (waitpid(-1, &stat, WNOHANG) > 0) {
	/* Subprocess exited; check status */
	if (!WIFEXITED(stat) \|\| WEXITSTATUS(stat)) {
	/* Exited with error; terminate siblings, clean up & bail out */
	for (proc = 0; proc < numProc; proc++)
	kill(procPids [proc], SIGKILL);

	/* Unmap shared memory, squish mapped file */
	munmap(photonCnt, sizeof(*photonCnt));
	close(shmFile);
	unlink(shmFname);
	error(USER, "failed photon distribution");
	}

	--proc;
	}

	/* Nod off for a bit and update progress */
	sleep(1);

	/* Asynchronous progress report from shared subprocess counters */
	shmLock(shmFile, F_RDLCK);
	repEmitted = repProgress = photonCnt -> numEmitted;
	repComplete = photonCnt -> numComplete;

	for (t = repProgress = repComplete = 0; t < NUM_PMAP_TYPES; t++)
	if ((pm = pmaps [t])) {
	/* Get global photon count from shmem updated by subprocs */
	repProgress += pm -> numPhotons = photonCnt -> numPhotons [t];
	repComplete += pm -> distribTarget;
	}

	repComplete *= numProc;
	shmLock(shmFile, F_UNLCK);

	if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
	pmapDistribReport();
	#ifdef SIGCONT
	else signal(SIGCONT, pmapDistribReport);
	#endif
	}
	#endif /* NIX */

	/* ===================================================================
	* POST-DISTRIBUTION - Set photon flux and build data struct for photon
	* storage, etc.
	* =================================================================== */
	#ifdef SIGCONT
	/* Reset signal handler */
	signal(SIGCONT, SIG_DFL);
	#endif
	free(emap.samples);

	/* Set photon flux */
	totalFlux /= photonCnt -> numEmitted;
	#if NIX
	/* Photon counters no longer needed, unmap shared memory */
	munmap(photonCnt, sizeof(*photonCnt));
	close(shmFile);
	unlink(shmFname);
	#else
	free(photonCnt);
	#endif
	if (verbose)
	eputs("\n");

	for (t = 0; t < NUM_PMAP_TYPES; t++)
	if (pmaps [t]) {
	if (verbose) {
	sprintf(errmsg, "Building %s photon map\n", pmapName [t]);
	eputs(errmsg);
	#if NIX
	fflush(stderr);
	#endif
	}

	/* Build underlying data structure; heap is destroyed */
	buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
	}

	/* Precompute photon irradiance if necessary */
	if (preCompPmap) {
	if (verbose)
	eputs("\n");
	preComputeGlobal(preCompPmap);
	}

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

pmap.cNo OneTemporaryActions

File Metadata

pmap.cView Options

Event Timeline

pmap.c
No OneTemporary
Actions

pmap.c
View Options