Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmapcontrib.c
a024a357a06cpmapmerge

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

RadiancePmap/pmapcontrib.ca024a357a06cpmapmerge

pmapcontrib.cView Options

RadiancePmap/pmapcontrib.c
a024a357a06cpmapmerge

pmapcontrib.c
View Options