pmapdata.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, May 23, 12:21

pmapdata.c
View Options

	#ifndef lint
	static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $";
	#endif

	/*
	=========================================================================
	Photon map types and high level interface to nearest neighbour lookups in
	underlying point cloud data structure.

	The default data structure is an in-core kd-tree (see pmapkdt.{h,c}).
	This can be overriden with the PMAP_OOC compiletime switch, which enables
	an out-of-core octree (see oococt.{h,c}).

	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",
	SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment")
	(c) Tokyo University of Science,
	supported by the JSPS Grants-in-Aid for Scientific Research
	(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
	=========================================================================

	$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $
	*/



	#include "pmapdata.h"
	#include "pmapdens.h"
	#include "pmaprand.h"
	#include "pmapmat.h"
	#include "pmaproi.h"
	#include "pmapcontrib.h"
	#include "otypes.h"
	#include "random.h"



	PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
	NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	#ifdef PMAP_PHOTONFLOW
	NULL, NULL,
	#endif
	NULL
	};



	/* Include routines to handle underlying point cloud data structure */
	#ifdef PMAP_OOC
	#include "pmapooc.c"
	#else
	#include "pmapkdt.c"
	#include "pmaptkdt.c"
	#endif

	/* Ambient include/exclude set (from ambient.c) */
	#ifndef MAXASET
	#define MAXASET 4095
	#endif
	extern OBJECT ambset [MAXASET+1];


	/* Callback to print photon attributes acc. to user defined format */
	int (printPhoton)(RAY r, Photon p, PhotonMap pm);



	/* PHOTON MAP BUILD ROUTINES ------------------------------------------ */

	void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
	/* Init photon map 'n' stuff... */
	{
	if (!pmap)
	return;

	pmap -> numPhotons = 0;
	pmap -> biasCompHist = NULL;
	pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE;
	pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE;
	pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0;
	pmap -> gatherTolerance = gatherTolerance;
	pmap -> minError = pmap -> maxError = pmap -> rmsError = 0;
	pmap -> numDensity = 0;
	pmap -> distribRatio = 1;
	pmap -> type = t;
	pmap -> squeue.node = NULL;
	pmap -> squeue.len = 0;
	#ifdef PMAP_PATHFILT
	/* Init path lookup table and its key buffer */
	pmap -> pathLUT = NULL;
	pmap -> pathLUTKeys = NULL;
	pmap -> numPathLUTKeys = 0;
	#endif

	/* Init stats */
	setcolor(pmap -> photonFlux, 0, 0, 0);
	pmap -> CoG [0] = pmap -> CoG [1] = pmap -> CoG [2] = 0;
	pmap -> CoGdist = 0;
	pmap -> maxDist2 = 0;

	/* Init transient photon mapping stuff */
	pmap -> velocity = photonVelocity;
	pmap -> minPathLen = pmap -> maxPathLen = pmap -> avgPathLen = 0;

	/* Init local RNG state */
	pmap -> randState [0] = 10243;
	pmap -> randState [1] = 39829;
	pmap -> randState [2] = 9433;
	pmapSeed(randSeed, pmap -> randState);

	/* Set up type-specific photon lookup callback */
	pmap -> lookup = pmapLookup [t];

	/* Init storage */
	pmap -> heap = NULL;
	pmap -> heapBuf = NULL;
	pmap -> heapBufLen = 0;

	/* Init precomputed contrib photon stuff */
	pmap -> contribMode = 0;
	pmap -> preCompContribTab = NULL;
	pmap -> contribHeap = NULL;
	pmap -> contribHeapBuf = NULL;
	pmap -> preCompContrib = NULL;
	pmap -> numContribSrc = 0;
	pmap -> contribSrc = NULL;
	/* Mark photon contribution source as unused */
	pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1;

	#ifdef PMAP_OOC
	OOC_Null(&pmap -> store);
	#else
	kdT_Null(&pmap -> store);
	#endif
	}



	void initPhotonHeap (PhotonMap *pmap)
	{
	int fdFlags;

	if (!pmap)
	error(INTERNAL, "undefined photon map in initPhotonHeap()");

	if (!pmap -> heap) {
	/* Open heap file */
	mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME));
	if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b")))
	error(SYSTEM, "failed opening heap file in initPhotonHeap()");

	#ifdef F_SETFL /* XXX is there an alternative needed for W33nd0z? */
	fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
	fcntl(fileno(pmap -> heap), F_SETFL, fdFlags \| O_APPEND);
	#endif /* ftruncate(fileno(pmap -> heap), 0); */
	}

	#ifdef PMAP_CONTRIB
	if (PMAP_CONTRIB_BINNING(pmap)) {
	/* Allocate precomputed contribution heap buffer */
	if (!pmap -> contribHeap) {
	/* Open heap file for precomputed contributions */
	mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME));
	if (!(pmap -> contribHeap = fopen(pmap -> contribHeapFname, "w+b")))
	error(SYSTEM, "failed opening precomputed contribution "
	"heap file in initPhotonHeap()"
	);
	#ifdef F_SETFL /* XXX is there an alternative needed for W33nd0z? */
	fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL);
	fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags \| O_APPEND);
	#endif /* ftruncate(fileno(pmap -> contribHeap), 0); */
	}
	}
	#endif
	}



	void flushPhotonHeap (PhotonMap *pmap)
	{
	int fd, fdContrib = -1;
	unsigned long len, lenContrib = 0;

	if (!pmap)
	error(INTERNAL, "undefined photon map in flushPhotonHeap()");

	if (!pmap -> heap \|\| !pmap -> heapBuf) {
	/* Silently ignore undefined heap
	error(INTERNAL, "undefined photon heap in flushPhotonHeap()"); */
	return;
	}

	fd = fileno(pmap -> heap);
	len = pmap -> heapBufLen * sizeof(Photon);

	#ifdef PMAP_CONTRIB
	if (PMAP_CONTRIB_BINNING(pmap)) {
	if (!pmap -> contribHeap \|\| !pmap -> contribHeapBuf) {
	/* Silently ignore undefined heap
	error(INTERNAL, "undefined contribution heap in flushPhotonHeap()");
	*/
	return;
	}

	/* Do same for the contribution photon buffa if binning is enabled.
	* Note the photon and contribution heap contain the same number of
	* entries to ensure photons are flushed in sync along with their
	* corresponding contributions to their respective heap files */
	fdContrib = fileno(pmap -> contribHeap);
	lenContrib = pmap -> heapBufLen *
	((PreComputedContrib*)pmap -> preCompContrib) -> contribSize;
	}
	#endif

	#if NIX
	/* Lock heap file(s) prior to writing */
	shmLock(fd, F_WRLCK);
	if (lenContrib)
	shmLock(fdContrib, F_WRLCK);
	#endif

	#ifdef DEBUG_PMAP
	sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n",
	getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon)
	);
	eputs(errmsg);
	#endif

	/* Atomically write write block to photon heap file */
	/* !!! Unbuffered I/O via pwrite() avoids potential race conditions
	* !!! and buffer corruption which can occur with lseek()/fseek()
	* !!! followed by write()/fwrite(). */
	/if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) /
	if (write(fd, pmap -> heapBuf, len) != len) {
	/* Clean up temp file */
	fclose(pmap -> heap);
	unlink(pmap -> heapFname);
	error(SYSTEM, "can't append to photon heap file in flushPhotonHeap()");
	}

	#ifdef PMAP_CONTRIB
	/* Atomically write block to contribution heap file */
	if (fdContrib > -1 &&
	write(fdContrib, pmap -> contribHeapBuf, lenContrib) != lenContrib
	) {
	/* Clean up temp file */
	fclose(pmap -> contribHeap);
	unlink(pmap -> contribHeapFname);
	error(SYSTEM, "can't append to contrib heap file in flushPhotonHeap()");
	}
	#endif

	#if NIX
	if (fsync(fd) \|\| fdContrib > -1 && fsync(fdContrib))
	error(SYSTEM, "failed to fsync heap in flushPhotonHeap()");

	/* Release lock on heap file(s) */
	shmLock(fd, F_UNLCK);
	if (fdContrib > -1)
	shmLock(fdContrib, F_UNLCK);
	#endif

	pmap -> heapBufLen = 0;
	}



	#ifdef DEBUG_PMAP
	static int checkPhotonHeap (FILE *file)
	/* Check heap for nonsensical or duplicate photons */
	{
	Photon p, lastp;
	int i, dup;

	rewind(file);
	memset(&lastp, 0, sizeof(lastp));

	while (fread(&p, sizeof(p), 1, file)) {
	dup = 1;

	for (i = 0; i <= 2; i++) {
	if (p.pos [i] < thescene.cuorg [i] \|\|
	p.pos [i] > thescene.cuorg [i] + thescene.cusize
	) {
	sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
	p.pos [0], p.pos [1], p.pos [2]
	);
	error(WARNING, errmsg);
	}

	dup &= p.pos [i] == lastp.pos [i];
	}

	if (dup) {
	sprintf(errmsg, "consecutive duplicate photon in heap "
	"at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2]
	);
	error(WARNING, errmsg);
	}
	}

	return 0;
	}
	#endif



	int newPhoton (PhotonMap* pmap, const RAY* ray, const char *contribs)
	{
	unsigned i;
	Photon photon;
	COLOR photonFlux;

	/* Account for distribution ratio */
	if (!pmap \|\| pmapRandom(pmap -> randState) > pmap -> distribRatio)
	return -1;

	/* Don't store on sources */
	if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
	return -1;

	/* Ignore photon if modifier in/outside exclude/include set */
	if (ambincl != -1 && ray -> ro &&
	ambincl != inset(ambset, ray->ro->omod)
	)
	return -1;

	/* Store photon if within a region of interest (for ze Ecksperts!) */
	if (!photonInROI(ray))
	return -1;

	/* Adjust flux according to distribution ratio and ray weight */
	copycolor(photonFlux, ray -> rcol);
	#if 0
	/* Factored out ray -> rweight as deprecated (?) for pmap, and infact
	erroneously attenuates volume photon flux based on extinction,
	which is already factored in by photonParticipate() */
	scalecolor(photonFlux,
	ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1)
	);
	#else
	scalecolor(photonFlux,
	1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1)
	);
	#endif
	setPhotonFlux(&photon, photonFlux);

	/* Set photon position and flags */
	VCOPY(photon.pos, ray -> rop);
	photon.flags = 0;
	photon.caustic = PMAP_CAUSTICRAY(ray);
	/* Set photon's subprocess index (to consolidate contribution sources
	after photon distribution completes) */
	photon.proc = PMAP_GETRAYPROC(ray);

	switch (pmap -> type) {
	#ifdef PMAP_CONTRIB
	case PMAP_TYPE_CONTRIB:
	/* Contrib photon before precomputation (i.e. in forward
	pass); set contribution source from last index in contrib
	source array. Note the contribution source bin has already
	been set by newPhotonContribSrc(). */
	photon.aux.contribSrc = pmap -> numContribSrc;
	/* Photon will be stored; set contribution source index,
	* thereby marking it as having spawned photon(s) */
	if (ray -> rsrc > PMAP_MAXSRCIDX)
	error(INTERNAL, "contribution source index overflow");
	else pmap -> lastContribSrc.srcIdx = ray -> rsrc;
	break;

	case PMAP_TYPE_CONTRIB_CHILD:
	/* Contrib photon being precomputed; set contribution
	source from index passed in ray.
	NOTE: This is currently only for info and consistency with the
	above, but not used by the lookup routines */
	if (ray -> rsrc > PMAP_MAXSRCIDX)
	error(INTERNAL, "contribution source index overflow");
	else photon.aux.contribSrc = ray -> rsrc;
	break;
	#endif
	case PMAP_TYPE_TRANSIENT:
	#ifdef PMAP_PHOTONFLOW
	case PMAP_TYPE_TRANSLIGHTFLOW:
	#endif
	/* Set cumulative path length for transient photon, corresponding
	to its time of flight. The cumulative path length is obtained
	relative to the maximum path length photonMaxPathDist.
	The last intersection distance is subtracted as this isn't
	factored in until photonRay() generates a new photon ray. */
	photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot;
	break;

	default:
	/* Set photon's path ID from ray serial num
	(supplemented by photon.proc below) */
	photon.aux.pathID = ray -> rno;
	}

	/* Set normal */
	for (i = 0; i <= 2; i++)
	photon.norm [i] = 127.0 * (isVolumePmap(pmap)
	#ifdef PMAP_PHOTONFLOW
	\|\| isLightFlowPmap(pmap)
	#endif
	? ray -> rdir [i] : ray -> ron [i]
	);

	if (!pmap -> heapBuf) {
	/* Lazily allocate photon heap buffa */
	#if NIX
	/* Randomise buffa size to temporally decorellate flushes in
	* multiprocessing mode */
	srandom(randSeed + getpid());
	pmap -> heapBufSize = PMAP_HEAPBUFSIZE * 0.5 * (1 + frandom());
	#else
	/* Randomisation disabled for single processes on WIN; also useful
	* for reproducability during debugging */
	pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
	#endif

	#ifdef PMAP_CONTRIB
	if (contribs && pmap -> preCompContrib)
	/* Adapt buffa size to accommodate precomputed contributions */
	pmap -> heapBufSize *= (float)sizeof(Photon) /
	((PreComputedContrib*)pmap -> preCompContrib) -> contribSize;
	#endif
	if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
	error(SYSTEM,
	"out of memory allocating heap buffer in newPhoton()"
	);
	pmap -> heapBufLen = 0;
	}

	/* Photon initialised; now append to heap buffa */
	memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));

	#ifdef PMAP_CONTRIB
	if (contribs && pmap -> preCompContrib) {
	const unsigned long contribSize =
	((PreComputedContrib*)pmap -> preCompContrib) -> contribSize;

	if (!pmap -> contribHeapBuf) {
	/* Lazily allocate contribution heap buffa of same size as photon
	* heap buffa */
	pmap -> contribHeapBuf = calloc(pmap -> heapBufSize, contribSize);
	if (!pmap -> contribHeapBuf)
	error(SYSTEM, "out of memory allocating precomputed "
	"contribution heap buffer in newPhoton()"
	);
	}

	/* Append precomputed contribs to heap buffa */
	memcpy(pmap -> contribHeapBuf + pmap -> heapBufLen * contribSize,
	contribs, contribSize
	);
	}
	#endif

	if (++pmap -> heapBufLen >= pmap -> heapBufSize)
	/* Photon (and contrib, if applicable) heap buffa(s) full,
	flush to heap file(s) */
	flushPhotonHeap(pmap);

	pmap -> numPhotons++;

	/* Print photon attributes */
	if (printPhoton)
	/* Non-const kludge */
	printPhoton((RAY*)ray, &photon, pmap);

	return 0;
	}



	void buildPhotonMap (PhotonMap pmap, double photonFlux,
	PhotonContribSourceIdx *contribSrcOfs, unsigned nproc
	)
	{
	unsigned long n, nCheck = 0;
	unsigned i;
	Photon *p;
	COLOR flux;
	char nuHeapFname [sizeof(PMAP_TMPFNAME)];
	FILE *nuHeap;
	/* Need double here to reduce summation errors */
	double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0,
	avgPathLen = 0;
	FVECT d;

	if (!pmap)
	error(INTERNAL, "undefined photon map in buildPhotonMap()");

	/* Get number of photons from heapfile size */
	if (fseek(pmap -> heap, 0, SEEK_END) < 0)
	error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()");
	pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);

	if (!pmap -> numPhotons)
	error(INTERNAL, "empty photon map in buildPhotonMap()");

	if (!pmap -> heap)
	error(INTERNAL, "no heap in buildPhotonMap()");

	#ifdef DEBUG_PMAP
	eputs("Checking photon heap consistency...\n");
	checkPhotonHeap(pmap -> heap);

	sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons);
	eputs(errmsg);
	#endif

	/* Allocate heap buffa */
	if (!pmap -> heapBuf) {
	pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
	pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon));
	if (!pmap -> heapBuf)
	error(SYSTEM,
	"failed to allocate postprocessed photon heap in buildPhotonMap()"
	);
	}

	/* We REALLY don't need yet another @%&*! heap just to hold the scaled
	* photons, but can't think of a quicker fix... */
	mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME));
	if (!(nuHeap = fopen(nuHeapFname, "w+b")))
	error(SYSTEM,
	"failed to open postprocessed photon heap in buildPhotonMap()"
	);

	rewind(pmap -> heap);

	#ifdef DEBUG_PMAP
	eputs("Postprocessing photons...\n");
	#endif

	while (!feof(pmap -> heap)) {
	#ifdef DEBUG_PMAP
	printf("Reading %lu at %lu... ",
	pmap -> heapBufSize, ftell(pmap->heap)
	);
	#endif
	pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
	pmap -> heapBufSize, pmap -> heap
	);
	#ifdef DEBUG_PMAP
	printf("Got %lu\n", pmap -> heapBufLen);
	#endif

	if (ferror(pmap -> heap))
	error(SYSTEM, "failed to read photon heap in buildPhotonMap()");

	for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
	/* Update min and max pos and set photon flux */
	for (i = 0; i <= 2; i++) {
	if (p -> pos [i] < pmap -> minPos [i])
	pmap -> minPos [i] = p -> pos [i];
	else if (p -> pos [i] > pmap -> maxPos [i])
	pmap -> maxPos [i] = p -> pos [i];

	/* Update centre of gravity with photon position */
	CoG [i] += p -> pos [i];
	}

	if (isContribPmap(pmap) && contribSrcOfs)
	/* Linearise contrib photon origin index from subprocess index
	* using the per-subprocess offsets in contribSrcOfs */
	p -> aux.contribSrc += contribSrcOfs [p -> proc];
	else if (isTransientPmap(pmap)
	#ifdef PMAP_PHOTONFLOW
	\|\| isTransLightFlowPmap(pmap)
	#endif
	) {
	/* Update min and max path length */
	if (p -> aux.pathLen < pmap -> minPathLen)
	pmap -> minPathLen = p -> aux.pathLen;
	else if (p -> aux.pathLen > pmap -> maxPathLen)
	pmap -> maxPathLen = p -> aux.pathLen;

	/* Update avg photon path length */
	avgPathLen += p -> aux.pathLen;
	}

	/* Scale photon's flux (hitherto normalised to 1 over RGB); in
	* case of a contrib photon map, this is done per light source,
	* and photonFlux is assumed to be an array */
	getPhotonFlux(p, flux);

	if (photonFlux) {
	scalecolor(flux,
	#if 1
	photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0]
	#else
	photonFlux [0]
	#endif
	);
	setPhotonFlux(p, flux);
	}

	/* Update average photon flux; need a double here */
	addcolor(avgFlux, flux);
	}

	/* Write modified photons to new heap */
	fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);

	if (ferror(nuHeap)) {
	error(SYSTEM,
	"failed postprocessing photon flux in buildPhotonMap()"
	);
	/* Clean up temp files */
	fclose(pmap -> heap);
	fclose(nuHeap);
	unlink(pmap -> heapFname);
	unlink(nuHeapFname);
	}

	nCheck += pmap -> heapBufLen;
	}

	/* !!! Flush pending writes to new heap !!! */
	fflush(nuHeap);

	#ifdef DEBUG_PMAP
	if (nCheck < pmap -> numPhotons)
	error(INTERNAL, "truncated photon heap in buildPhotonMap");
	#endif

	/* Finalise average photon flux */
	scalecolor(avgFlux, 1.0 / pmap -> numPhotons);
	copycolor(pmap -> photonFlux, avgFlux);

	/* Average photon positions to get centre of gravity */
	for (i = 0; i < 3; i++)
	pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons;

	/* Average photon path lengths */
	pmap -> avgPathLen = avgPathLen /= pmap -> numPhotons;

	rewind(pmap -> heap);
	/* Compute average photon distance to centre of gravity */
	while (!feof(pmap -> heap)) {
	pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
	pmap -> heapBufSize, pmap -> heap
	);

	for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
	VSUB(d, p -> pos, CoG);
	CoGdist += DOT(d, d);

	if (isTransientPmap(pmap)
	#ifdef PMAP_PHOTONFLOW
	\|\| isTransLightFlowPmap(pmap)
	#endif
	) {
	/* Add temporal distance (in units of photon path length)
	for transient photons */
	d [0] = p -> aux.pathLen - avgPathLen;
	CoGdist += d [0] * d [0];
	}
	}
	}

	pmap -> CoGdist = CoGdist /= pmap -> numPhotons;

	/* Swap heaps, discarding unscaled photons */
	fclose(pmap -> heap);
	unlink(pmap -> heapFname);
	pmap -> heap = nuHeap;
	strcpy(pmap -> heapFname, nuHeapFname);

	#ifdef PMAP_OOC
	OOC_BuildPhotonMap(pmap, nproc);
	#else
	if (isTransientPmap(pmap)
	#ifdef PMAP_PHOTONFLOW
	\|\| isTransLightFlowPmap(pmap)
	#endif
	)
	kdT_TransBuildPhotonMap(pmap);
	else
	kdT_BuildPhotonMap(pmap);
	#endif

	/* Trash heap and its buffa */
	free(pmap -> heapBuf);
	fclose(pmap -> heap);
	unlink(pmap -> heapFname);
	pmap -> heap = NULL;
	pmap -> heapBuf = NULL;
	}



	/* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */

	/* Dynamic max photon search radius increase and reduction factors */
	#define PMAP_MAXDIST_INC 4
	#define PMAP_MAXDIST_DEC 0.9

	/* Num successful lookups before reducing in max search radius */
	#define PMAP_MAXDIST_CNT 1000

	/* Threshold below which we assume increasing max radius won't help */
	#define PMAP_SHORT_LOOKUP_THRESH 1

	/* Coefficient for adaptive maximum search radius */
	#define PMAP_MAXDIST_COEFF 100

	void findPhotons (PhotonMap* pmap, const RAY* ray)
	{
	int redo = 0;
	const RREAL norm = ray -> ron, pos = ray -> rop;

	if (!pmap -> squeue.len) {
	/* Lazy init priority queue */
	#ifdef PMAP_OOC
	OOC_InitFindPhotons(pmap);
	#else
	kdT_InitFindPhotons(pmap);
	#endif
	pmap -> minGathered = pmap -> maxGather;
	pmap -> maxGathered = pmap -> minGather;
	pmap -> totalGathered = 0;
	pmap -> numLookups = pmap -> numShortLookups = 0;
	pmap -> shortLookupPct = 0;
	pmap -> minError = FHUGE;
	pmap -> maxError = -FHUGE;
	pmap -> rmsError = 0;
	/* SQUARED max search radius limit is based on avg photon distance to
	* centre of gravity, unless fixed by user (maxDistFix > 0) */
	pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0
	? maxDistFix * maxDistFix
	: PMAP_MAXDIST_COEFF * pmap -> squeue.len *
	pmap -> CoGdist / pmap -> numPhotons;
	}

	do {
	pmap -> squeue.tail = 0;
	pmap -> maxDist2 = pmap -> maxDist0;

	if (isVolumePmap(pmap)
	#ifdef PMAP_PHOTONFLOW
	\|\| isLightFlowPmap(pmap) && sphericalIrrad
	#endif
	) {
	/* Volume photons are not associated with a normal or intersection
	point; null the normal and set origin as lookup point. */
	/* If a lightflow photon map is used and sphericalIrrad = 1,
	the mean spherical irradiance is evaluated, so normals are
	irrelevant and not filtered. */
	norm = NULL;
	pos = ray -> rorg;
	}

	/* XXX/NOTE: status returned by XXX_FindPhotons() is currently
	ignored; if no photons are found, an empty queue is returned
	under the assumption all photons are too distant to contribute
	significant flux. */
	#ifdef PMAP_OOC
	OOC_FindPhotons(pmap, pos, norm);
	#else
	if (isTransientPmap(pmap)
	#ifdef PMAP_PHOTONFLOW
	\|\| isTransLightFlowPmap(pmap)
	#endif
	)
	kdT_TransFindPhotons(pmap, pos, norm);
	else
	kdT_FindPhotons(pmap, pos, norm);
	#endif

	#ifdef PMAP_LOOKUP_INFO
	fprintf(stderr, "%d/%d %s photons found within radius %.3f "
	"at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail,
	pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2),
	ray -> rop [0], ray -> rop [1], ray -> rop [2],
	ray -> ro ? ray -> ro -> oname : "<null>"
	);
	#endif

	if (pmap -> squeue.tail < pmap -> squeue.len * pmap->gatherTolerance) {
	/* Short lookup; too few photons found */
	if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) {
	/* Ignore short lookups which return fewer than
	* PMAP_SHORT_LOOKUP_THRESH photons under the assumption there
	* really are no photons in the vicinity, and increasing the max
	* search radius therefore won't help */
	#ifdef PMAP_LOOKUP_WARN
	sprintf(errmsg,
	"%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
	pmap -> squeue.tail, pmap->squeue.len, pmapName [pmap->type],
	ray -> rop [0], ray -> rop [1], ray -> rop [2],
	ray -> ro ? ray -> ro -> oname : "<null>"
	);
	error(WARNING, errmsg);
	#endif

	/* Bail out after warning if maxDist is fixed */
	if (maxDistFix > 0)
	return;

	if (pmap -> maxDist0 < pmap -> maxDist2Limit) {
	/* Increase max search radius if below limit & redo search */
	pmap -> maxDist0 *= PMAP_MAXDIST_INC;
	#ifdef PMAP_LOOKUP_REDO
	redo = 1;
	#endif
	#ifdef PMAP_LOOKUP_WARN
	sprintf(errmsg, redo
	? "restarting photon lookup with max radius %.1e"
	: "max photon lookup radius adjusted to %.1e",
	pmap -> maxDist0
	);
	error(WARNING, errmsg);
	#endif
	}
	#ifdef PMAP_LOOKUP_REDO
	else {
	sprintf(errmsg, "max photon lookup radius clamped to %.1e",
	pmap -> maxDist0
	);
	error(WARNING, errmsg);
	}
	#endif
	}

	/* Reset successful lookup counter */
	pmap -> numLookups = 0;
	}
	else {
	/* Skip search radius reduction if maxDist is fixed */
	if (maxDistFix > 0)
	return;

	/* Increment successful lookup counter and reduce max search radius if
	* wraparound */
	pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT;
	if (!pmap -> numLookups)
	pmap -> maxDist0 *= PMAP_MAXDIST_DEC;

	redo = 0;
	}

	} while (redo);
	}



	Photon find1Photon (PhotonMap pmap, const RAY* ray, Photon *photon,
	void *photonIdx
	)
	{
	/* Init (squared) search radius to avg photon dist to centre of gravity */
	float maxDist2_0 = pmap -> CoGdist;
	int found = 0;
	#ifdef PMAP_LOOKUP_REDO
	#define REDO 1
	#else
	#define REDO 0
	#endif

	do {
	pmap -> maxDist2 = maxDist2_0;
	#ifdef PMAP_OOC
	found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon,
	(OOC_DataIdx*)photonIdx
	);
	#else
	found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
	#endif
	if (found < 0) {
	/* Expand search radius to retry */
	maxDist2_0 *= 2;
	#ifdef PMAP_LOOKUP_WARN
	sprintf(errmsg, "failed 1-NN photon lookup"
	#ifdef PMAP_LOOKUP_REDO
	", retrying with search radius %.2f", maxDist2_0
	#endif
	);
	error(WARNING, errmsg);
	#endif
	}
	} while (REDO && found < 0);

	/* Return photon buffer containing valid photon, else NULL */
	return found < 0 ? NULL : photon;
	}



	void getPhoton (PhotonMap pmap, PhotonIdx idx, Photon photon)
	{
	#ifdef PMAP_OOC
	if (OOC_GetPhoton(pmap, idx, photon))
	#else
	if (kdT_GetPhoton(pmap, idx, photon))
	#endif
	error(INTERNAL, "failed photon lookup");
	}



	Photon getNearestPhoton (const PhotonSearchQueue squeue, PhotonIdx idx)
	{
	#ifdef PMAP_OOC
	return OOC_GetNearestPhoton(squeue, idx);
	#else
	return kdT_GetNearestPhoton(squeue, idx);
	#endif
	}



	PhotonIdx firstPhoton (const PhotonMap *pmap)
	{
	#ifdef PMAP_OOC
	return OOC_FirstPhoton(pmap);
	#else
	return kdT_FirstPhoton(pmap);
	#endif
	}



	/* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */

	void deletePhotons (PhotonMap* pmap)
	{
	#ifdef PMAP_OOC
	OOC_Delete(&pmap -> store);
	#else
	kdT_Delete(&pmap -> store);
	#endif

	free(pmap -> squeue.node);
	free(pmap -> biasCompHist);

	pmap -> numPhotons = pmap -> minGather = pmap -> maxGather =
	pmap -> squeue.len = pmap -> squeue.tail = 0;

	#ifdef PMAP_PATHFILT
	if (pmap -> pathLUT) {
	lu_done(pmap -> pathLUT);
	free(pmap -> pathLUT);
	}
	if (pmap -> pathLUTKeys) {
	unsigned k;

	for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++]));
	free(pmap -> pathLUTKeys);
	}
	#endif

	/* Contrib stuff */
	if (isContribPmap(pmap) && pmap -> preCompContribTab) {
	/* Parent contrib photon map; clean up child photon maps containing
	precomputed contribs via lu_done() */
	lu_done(pmap -> preCompContribTab);
	free(pmap -> preCompContribTab);
	pmap -> preCompContribTab = NULL;

	if (pmap -> contribHeapBuf) {
	free(pmap-> contribHeapBuf);
	pmap -> contribHeapBuf = NULL;
	}
	if (pmap -> contribHeap) {
	fclose(pmap -> contribHeap);
	pmap -> contribHeap = NULL;
	}
	}
	}

pmapdata.cNo OneTemporaryActions

File Metadata

pmapdata.cView Options

Event Timeline

pmapdata.c
No OneTemporary
Actions

pmapdata.c
View Options