pmapcontrib.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Apr 27, 13:36

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 handle contribution binning, compression and encoding,
	and are used by 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 "pmapdiag.h"
	#include "pmaprand.h"
	#include "pmapmat.h"
	#include "pmapsrc.h"
	#include "pmutil.h"
	#include "otspecial.h"
	#include "otypes.h"
	#include "lookup.h"



	#ifdef PMAP_CONTRIB
	/* The following are convenient placeholders interfacing to mkpmap
	and rcontrib. They can be externally set via initPmapContribTab()
	and referenced within the contrib pmap modules. These variables
	can then be safely ignored by rtrace/rpict/rvu, without annoying
	linking errors. */
	/* Global pointer to rcontrib's contribution binning LUT */
	LUTAB *pmapContribTab = NULL;
	/* Contribution/coefficient mode flag */
	int *pmapContribMode;



	extern void SDdisk2square(double sq[2], double diskx, double disky);



	static int ray2bin (const RAY *ray, unsigned scDim)
	/* Map ray dir (pointing away from origin) to its 1D bin in an
	(scDim x scDim) Shirley-Chiu square.
	Returns -1 if mapped bin is invalid (e.g. behind plane) */
	{
	static int scRHS, varInit = 0;
	static FVECT scNorm, scUp;
	unsigned scBin [2];
	FVECT diskPlane;
	RREAL dz, diskx, disky, rad, diskd2, scCoord [2];

	if (!varInit) {
	/* Lazy init shirley-Chiu mapping orientation from function variables */
	scRHS = varvalue(PMAP_CONTRIB_SCRHS);
	scNorm [0] = varvalue(PMAP_CONTRIB_SCNORMX);
	scNorm [1] = varvalue(PMAP_CONTRIB_SCNORMY);
	scNorm [2] = varvalue(PMAP_CONTRIB_SCNORMZ);
	scUp [0] = varvalue(PMAP_CONTRIB_SCUPX);
	scUp [1] = varvalue(PMAP_CONTRIB_SCUPY);
	scUp [2] = varvalue(PMAP_CONTRIB_SCUPZ);
	varInit ^= 1;
	}

	/* Map incident direction to disk */
	dz = DOT(ray -> rdir, scNorm); /* normal proj */
	if (dz > 0) {
	fcross(diskPlane, scUp, scNorm); /* y-axis in plane, perp to up */
	diskPlane [0] *= scRHS;
	diskx = DOT(ray -> rdir, diskPlane);
	disky = DOT(ray -> rdir, scUp) - dz * DOT(ray -> rdir, scUp);
	diskd2 = diskx * diskx + disky * disky; /* in-plane length^2 */
	rad = dz>FTINY ? sqrt((1 - dzdz) / diskd2) : 0; / radial factor */
	disky = rad; / diskx, disky now in range [-1, 1] */

	/* Apply Shirley-Chiu mapping of (diskx, disky) to square */
	SDdisk2square(scCoord, diskx, disky);
	/* Map Shirley-Chiu square coords to 1D bin */
	scBin [0] = scCoord [0] * scDim;
	scBin [1] = scCoord [1] * scDim;
	return PMAP_CONTRIB_XY2LIN(scDim, scBin [0], scBin [1]);
	}
	else return -1;
	}


	/* ------------------ CONTRIBSRC STUFF --------------------- */


	#ifndef PMAP_CONTRIB_TEST
	MODCONT addContribModifier (LUTAB contribTab, unsigned *numContribs,
	char mod, char binParm, int binCnt
	)
	/* Add light source modifier mod to contribution lookup table contribTab,
	and update numContribs. Return initialised contribution data for this
	modifier. This code adapted from rcontrib.c:addmodifier(). */
	{
	LUENT *lutEntry = lu_find(contribTab, mod);
	MODCONT *contrib;
	EPNODE *eBinVal;
	unsigned numCoeffs;

	/* Warn if potential coefficient index overflow in mRGBE encoding.
	This requires getting the number of wavelet coefficients generated
	by the transform a priori. */
	if (!(numCoeffs = padWaveletXform2(NULL, NULL, sqrt(binCnt), NULL))) {
	sprintf(errmsg, "can't determine number of wavelet coefficients "
	"for modifier %s", mod
	);
	error(INTERNAL, errmsg);
	}

	if (numCoeffs * numCoeffs > PMAP_CONTRIB_MAXCOEFFS) {
	sprintf(errmsg, "mRGBE data field may overflow for modifier %s; "
	"reduce -bn if precomputing contributions fails", mod
	);
	error(WARNING, errmsg);
	}

	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 (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, 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, binCnt
	);
	}



	static int contribSourceBin (LUTAB contribs, const RAY ray)
	/* Map contribution source ray to its bin for light source ray -> rsrc,
	using Shirley-Chiu disk-to-square mapping.
	Return invalid bin -1 if mapping failed. */
	{
	const SRCREC *src;
	const OBJREC *srcMod;
	const MODCONT *srcCont;
	RAY srcRay;
	int bin, i;

	/* 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)
	/* Not interested in this source (modifier not in contrib LUT) */
	return -1;

	/* Set up shadow ray pointing to source for disk2square mapping */
	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;
	#if 0
	worldfunc(RCCONTEXT, &srcRay);
	#endif
	set_eparams((char*)srcCont -> params);

	if ((bin = ray2bin(&srcRay, sqrt(srcCont -> nbins))) < 0)
	error(WARNING, "Ignoring invalid bin in contribSourceBin()");

	return bin;
	}



	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 contrib
	* 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
	)
	error(INTERNAL, "contribution source overflow");
	}

	if (contribSrcRay) {
	/* 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 contrib source will be
	used, since the ray will be lost once a photon is spawned */
	pmap -> lastContribSrc.srcBin = contribSourceBin(
	pmapContribTab, contribSrcRay
	);

	if (pmap -> lastContribSrc.srcBin < 0) {
	/* Warn if invalid bin, but trace photon nonetheless. It will
	count as emitted to prevent bias, but will not be stored in
	newPhoton(), as it contributes zero flux */
	sprintf(errmsg, "invalid bin for light source %s, "
	"contribution photons will be discarded",
	source [contribSrcRay -> rsrc].so -> oname
	);
	error(WARNING, errmsg);
	}
	}

	return pmap -> numContribSrc;
	}



	PhotonContribSourceIdx buildContribSources (PhotonMap *pmap,
	FILE contribSrcHeap, char contribSrcHeapFname,
	PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps
	)
	/* Consolidate per-subprocess contribution source 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 contrib 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;
	}



	/* ----------------- CONTRIB PMAP ENCODING STUFF -------------------- */

	static void coeffSwap (PreComputedContribCoeff *c1,
	PreComputedContribCoeff *c2
	) {
	PreComputedContribCoeff tCoeff;

	tCoeff.coeff = c1 -> coeff;
	tCoeff.idx = c1 -> idx;
	c1 -> coeff = c2 -> coeff;
	c1 -> idx = c2 -> idx;
	c2 -> coeff = tCoeff.coeff;
	c2 -> idx = tCoeff.idx;
	}



	static int coeffPartition (PreComputedContribCoeff *coeffs,
	unsigned medianPos, unsigned left, unsigned right
	)
	/* REVERSE partition coefficients by magnitude, such that
	coeffs [left..medianPos] >= coeffs [medianPos+1..right].
	Returns median's position. */
	{
	unsigned long l, r, m;
	WAVELET_COEFF lVal, rVal, mVal, tVal;

	#define COEFFVAL(c) (DOT((c) -> coeff, (c) -> coeff))

	if (left < right) {
	/* Select pivot from median-of-three and move to photons [right]
	(a.k.a. Lomuto partitioning) */
	l = left;
	r = right;
	m = l + ((r - l) >> 1); /* Avoids overflow vs. (l+r) >> 1 */
	lVal = COEFFVAL(coeffs + l);
	rVal = COEFFVAL(coeffs + r);
	mVal = COEFFVAL(coeffs + m);
	if (mVal > lVal) {
	coeffSwap(coeffs + m, coeffs + l);
	tVal = mVal;
	mVal = lVal;
	lVal = tVal;
	}
	if (rVal > lVal) {
	coeffSwap(coeffs + r, coeffs + l);
	tVal = rVal;
	rVal = lVal;
	lVal = tVal;
	}
	if (mVal > rVal) {
	coeffSwap(coeffs + m, coeffs + r);
	tVal = mVal;
	mVal = rVal;
	rVal = tVal;
	}
	/* Pivot with key rVal is now in coeffs [right] */

	/* l & r converge, swapping coefficients out of (reversed) order
	with respect to pivot. The convergence point l = r is the
	pivot's final position */
	while (l < r) {
	while (l < r && COEFFVAL(coeffs + l) > rVal)
	l++;
	while (r > l && COEFFVAL(coeffs + r) <= rVal)
	r--;

	/* Coeffs out of order, must swap */
	if (l < r)
	coeffSwap(coeffs + l, coeffs + r);
	};

	/* Now l == r is pivot's final position; swap these coeffs */
	coeffSwap(coeffs + l, coeffs + right);

	/* Recurse in partition containing median */
	if (l > medianPos)
	return coeffPartition(coeffs, medianPos, left, l - 1);
	else if (l < medianPos)
	return coeffPartition(coeffs, medianPos, l + 1, right);
	else /* l == medianPos, partitioning done */
	return l;
	}
	else return left;
	}



	static int coeffIdxCompare (const void c1, const void c2)
	/* Comparison function to sort thresholded coefficients by index */
	{
	const PreComputedContribCoeff tcoeff1 = c1, tcoeff2 = c2;
	const unsigned v1 = tcoeff1 -> idx,
	v2 = tcoeff2 -> idx;

	if (v1 < v2)
	return -1;
	else if (v1 > v2)
	return 1;
	else
	return 0;
	}



	static int thresholdContribs (PreComputedContrib *preCompContrib)
	/* Threshold wavelet detail coefficients in preCompContrib ->
	waveletMatrix [approxDim..coeffDim-1] [approxDim..coeffDim-1] (where
	approxDim = WAVELET_PADD4_APPROXDIM) by keeping the (preCompContrib ->
	nCompressedCoeffs) largest of these and returning them in
	preCompContrib -> compressedCoeffs along with their original matrix
	indices.
	NOTE: The wavelet approximation coefficients
	preCompContrib -> waveletMatrix [0..approxDim-1] [0..approxDim-1]
	and excluded from thresholding to minimise compression artefacts.
	Returns 0 on success, else -1. */
	{
	unsigned i, j, coeffDim, coeffIdx, nNzDetailCoeffs,
	nCompressedCoeffs, numThresh;
	WaveletMatrix2 waveletMatrix;
	PreComputedContribCoeff threshCoeffs, threshCoeffPtr;

	if (!preCompContrib \|\| !(coeffDim = preCompContrib -> coeffDim) \|\|
	!(threshCoeffs = preCompContrib -> compressedCoeffs) \|\|
	!(waveletMatrix = preCompContrib -> waveletMatrix)
	)
	/* Should be initialised by caller! */
	return -1;

	/* Set up coefficient buffer for compression (thresholding), skipping
	* the approximation coefficients in the upper left of waveletMatrix,
	* which are not thresholded. Also skip zero coefficients (resulting
	* from padding), since these are already implicitly thresholded. The
	* original 2D matrix indices are linearised to 1D and saved with each
	* coefficient to restore the original sparse coefficient matrix. */
	for (i = 0, threshCoeffPtr = threshCoeffs; i < coeffDim; i++)
	for (j = 0; j < coeffDim; j++)
	if ((i >= WAVELET_PADD4_APPROXDIM \|\|
	j >= WAVELET_PADD4_APPROXDIM
	) && !coeffIsZero(waveletMatrix [i] [j])
	) {
	/* Nonzero detail coefficient; set up pointer to coeff
	(sorts faster than 3 doubles) and save original
	(linearised) matrix index */
	threshCoeffPtr -> idx = PMAP_CONTRIB_XY2LIN(coeffDim, i, j);
	threshCoeffPtr++ -> coeff = (WAVELET_COEFF*)&(
	waveletMatrix [i] [j]
	);
	}

	/* Num of nonzero detail coeffs in buffer, number actually expected */
	numThresh = threshCoeffPtr - threshCoeffs;
	nNzDetailCoeffs = WAVELET_PADD4_NUMDETAIL(
	preCompContrib -> nNonZeroCoeffs
	);
	nCompressedCoeffs = preCompContrib -> nCompressedCoeffs;

	/* If fewer nonzero detail coeff are in the threshold buffer than
	* anticipated, the loop below fills the remainder of the threshold
	* buffer with duplicates of a coefficient in the lower right of the
	* matrix, which is padding and guaranteed to be zero. This condition
	* can occur if the wavelet transform actually generated genuine zero
	* detail coefficients. Infact it's quite common if the wavelet
	* transformed contributions are quite sparse. */
	for (i = j = coeffDim - 1; numThresh < nNzDetailCoeffs; numThresh++) {
	threshCoeffPtr -> idx = PMAP_CONTRIB_XY2LIN(coeffDim, i, j);
	threshCoeffPtr++ -> coeff = (WAVELET_COEFF*)(waveletMatrix [i][j]);
	}

	/* Partition coeffs in reverse order, such that all coeffs in
	threshCoeffs [0..nCompressedCoeffs-1] are larger than those in
	threshCoeffs [nCompressedCoeffs..nNzDetailCoeffs-1] */
	coeffPartition(threshCoeffs, nCompressedCoeffs - 1,
	0, nNzDetailCoeffs - 1
	);

	#ifdef PMAP_CONTRIB_DBG
	/* Check coefficient partitioning */
	threshCoeffPtr = threshCoeffs + nCompressedCoeffs - 1;
	for (i = 0; i < nCompressedCoeffs; i++)
	if (COEFFVAL(threshCoeffs + i) < COEFFVAL(threshCoeffPtr))
	error(CONSISTENCY, "invalid wavelet coefficient partitioning "
	"in thresholdContribs()"
	);

	for (; i < nNzDetailCoeffs; i++)
	if (COEFFVAL(threshCoeffs + i) > COEFFVAL(threshCoeffPtr))
	error(CONSISTENCY, "invalid wavelet coefficient partitioning "
	"in thresholdContribs()"
	);
	#endif

	/* Sort partition containing nCompressedCoeffs coefficients by index
	* (ignoring the remaining coefficients since these are now dropped
	* due to tresholding) */
	qsort(threshCoeffs, nCompressedCoeffs,
	sizeof(PreComputedContribCoeff), coeffIdxCompare
	);

	return 0;
	}



	static int encodeContribs (PreComputedContrib *preCompContrib,
	float compressRatio
	)
	/* Apply wavelet transform to input matrix preCompContrib ->
	waveletMatrix and compress according to compressRatio, storing
	thresholded and mRGBE-encoded coefficients in preCompContrib ->
	mrgbeCoeffs. Note that the approximation coefficients in the upper
	left of the matrix are not encoded, and returned in
	preCompContrib -> waveletMatrix
	[0..WAVELET_PADD4_APPROXDIM-1] [0..WAVELET_PADD4_APPROXDIM-1].
	Returns 0 on success. */
	{
	unsigned i, j, k, scDim, lastCoeffIdx;
	WaveletMatrix2 waveletMatrix, tWaveletMatrix;
	PreComputedContribCoeff *threshCoeffs;
	mRGBERange *mrgbeRange;
	mRGBE *mrgbeCoeffs;
	WAVELET_COEFF absCoeff;
	#ifdef PMAP_CONTRIB_DBG
	WaveletCoeff3 decCoeff;
	unsigned decIdx;
	#endif

	if (!preCompContrib \|\| !preCompContrib -> mrgbeCoeffs \|\|
	!preCompContrib -> compressedCoeffs \|\|
	!(waveletMatrix = preCompContrib -> waveletMatrix) \|\|
	!(tWaveletMatrix = preCompContrib -> tWaveletMatrix) \|\|
	!(scDim = preCompContrib -> scDim)
	)
	/* Should be initialised by the caller! */
	return -1;

	#ifdef PMAP_CONTRIB_DBG
	for (i = 0; i < scDim; i++) {
	for (j = 0; j < scDim; j++) {
	for (k = 0; k < 3; k++) {
	#if 0
	/* Set contributions to bins for debugging */
	waveletMatrix [i] [j] [k] = PMAP_CONTRIB_XY2LIN(scDim, i, j);
	#elif 0
	/* Replace contribs with "bump" function */
	waveletMatrix [i] [j] [k] =
	(1. - fabs(i - scDim/2. + 0.5) / (scDim/2. - 0.5)) *
	(1. - fabs(j - scDim/2. + 0.5) / (scDim/2. - 0.5));
	#endif
	}

	#if 0
	/* Dump contribs prior to encoding for debugging */
	printf("% 7.3g\t", colorAvg(waveletMatrix [i] [j]));
	}
	putchar('\n');
	}
	putchar('\n');
	#else
	}
	}
	#endif
	#endif

	/* Do 2D wavelet transform on preCompContrib -> waveletMatrix */
	if (padWaveletXform2(waveletMatrix, tWaveletMatrix, scDim, NULL) !=
	preCompContrib -> coeffDim
	)
	error(INTERNAL, "failed wavelet transform in encodeContribs()");

	/* Compress wavelet detail coeffs by thresholding */
	if (thresholdContribs(preCompContrib) < 0)
	error(INTERNAL, "failed wavelet compression in encodeContribs()");
	threshCoeffs = preCompContrib -> compressedCoeffs;

	/* Init per-channel coefficient range for mRGBE encoding */
	mrgbeRange = &preCompContrib -> mrgbeRange;
	setcolor(mrgbeRange -> min, FHUGE, FHUGE, FHUGE);
	setcolor(mrgbeRange -> max, 0, 0, 0);

	/* Update per-channel coefficient range */
	for (i = 0; i < preCompContrib -> nCompressedCoeffs; i++) {
	for (k = 0; k < 3; k++) {
	#ifdef PMAP_CONTRIB_DBG
	#if 0
	/* Replace wavelet coeff with its linear index for debugging */
	threshCoeffs [i].coeff [k] = threshCoeffs [i].idx = i + 1;
	#endif
	#endif
	absCoeff = fabs(threshCoeffs [i].coeff [k]);

	if (absCoeff < mrgbeRange -> min [k])
	mrgbeRange -> min [k] = absCoeff;
	if (absCoeff > mrgbeRange -> max [k])
	mrgbeRange -> max [k] = absCoeff;
	}
	}

	if (preCompContrib -> nCompressedCoeffs == 1)
	/* Maximum compression with just 1 (!) compressed detail coeff (is
	* this even useful?), so mRGBE range is undefined since min & max
	* coincide; set minimum to 0, maximum to the single remaining
	* coefficient */
	setcolor(mrgbeRange -> min, 0, 0, 0);

	/* Init mRGBE coefficient normalisation from range */
	if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max))
	return -1;

	mrgbeCoeffs = preCompContrib -> mrgbeCoeffs;

	/* Encode wavelet detail coefficients to mRGBE */
	for (i = lastCoeffIdx = 0; i < preCompContrib -> nCompressedCoeffs;
	lastCoeffIdx = threshCoeffs [i++].idx
	) {
	/* HACK: To reduce the likelihood of overflowing the mRGBE data
	* field with the coefficient index, it is expressed incrementally
	* w.r.t. the previously encoded coefficient's index, instead of
	* absolutely. This implies threshCoeffs must be sorted by
	* coefficient index to ensure increments are positive, and to
	* minimise their magnitude.
	* Note that an overflow cannot be predicted beforehand, e.g. by
	* mkpmap when parsing the number of bins, as this depends on the
	* sparseness of the wavelet coefficients (which in turn depends on
	* the frequency distribution of the binned contributions), and the
	* fraction of those that are dropped (i.e. the compression ratio).
	* */
	mrgbeCoeffs [i] = mRGBEencode(threshCoeffs [i].coeff,
	mrgbeRange, threshCoeffs [i].idx - lastCoeffIdx
	);
	if (!mrgbeCoeffs [i].all)
	error(INTERNAL, "failed mRGBE encoding in encodeContribs()");
	#ifdef PMAP_CONTRIB_DBG
	/* Encoding sanity check */
	decIdx = mRGBEdecode(mrgbeCoeffs [i], mrgbeRange, decCoeff);
	#if 1
	if (decIdx != threshCoeffs [i].idx - lastCoeffIdx \|\|
	sqrt(dist2(decCoeff, threshCoeffs [i].coeff)) >
	0.1 * colorAvg(mrgbeRange -> max)
	) {
	sprintf(errmsg, "failed sanity check in encodeContribs()\n"
	"Encoded: [%.3g %.3g %.3g %d]\nDecoded: [%.3g %.3g %.3g %d]",
	threshCoeffs [i].coeff [0], threshCoeffs [i].coeff [1],
	threshCoeffs [i].coeff [2], threshCoeffs [i].idx,
	decCoeff [0], decCoeff [1], decCoeff [2],
	decIdx - lastCoeffIdx
	);
	error(CONSISTENCY, errmsg);
	}
	#endif
	for (k = 0; k < 3; k++)
	if (decCoeff [k] / threshCoeffs [i].coeff [k] < 0)
	error(CONSISTENCY,
	"coefficient sign reversal in encodeContribs()"
	);
	#endif
	}

	return 0;
	}



	static void initContribHeap (PhotonMap *pmap)
	/* Initialise precomputed contribution heap */
	{
	int fdFlags;

	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 initContribHeap()"
	);
	#ifdef F_SETFL /* XXX is there an alternative needed for Wind0z? */
	fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL);
	fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags \| O_APPEND);
	#endif /* ftruncate(fileno(pmap -> heap), 0); */
	}
	}



	static MODCONT getPhotonContrib (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(
	pmapContribTab, srcMod -> oname
	) -> data;
	Photon *photon;
	COLOR flux;
	PhotonSearchQueueNode *sqn;
	float r2, norm;
	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;
	}

	/* Avg radius^2 between furthest two photons to improve accuracy */
	sqn = pmap -> squeue.node + 1;
	r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
	r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));

	/* XXX: OMIT extra normalisation factor 1/PI for ambient calc? */
	#ifdef PMAP_EPANECHNIKOV
	/* Normalise accumulated flux by Epanechnikov kernel integral in 2D
	(see Eq. 4.4, p.76 in Silverman, "Density Estimation for
	Statistics and Data Analysis", 1st Ed., 1986, and
	Wann Jensen, "Realistic Image Synthesis using Photon Mapping"),
	include RADIANCE-specific lambertian factor PI */
	norm = 2 / (PI * PI * r2);
	#else
	/* Normalise accumulated flux by search area PI * r^2, including
	RADIANCE-specific lambertian factor PI */
	norm = 1 / (PI * PI * r2);
	#endif

	/* 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, norm);
	#ifdef PMAP_EPANECHNIKOV
	/* Apply Epanechnikov kernel to photon flux based on distance */
	scalecolor(flux, 1 - sqn -> dist2 / r2);
	#endif
	addcolor(irrad, flux);
	addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux);
	}

	return contrib;
	}



	void freePreCompContribNode (void *p)
	/* Clean up precomputed contribution LUT entry */
	{
	PhotonMap preCompContribPmap = (PhotonMap)p;
	PreComputedContrib *preCompContrib;

	if (preCompContribPmap) {
	preCompContrib = (PreComputedContrib*)(
	preCompContribPmap -> preCompContrib
	);

	if (preCompContrib) {
	/* Free primary and transposed wavelet matrices */
	freeWaveletMatrix2(preCompContrib -> waveletMatrix,
	preCompContrib -> coeffDim
	);
	freeWaveletMatrix2(preCompContrib -> tWaveletMatrix,
	preCompContrib -> coeffDim
	);

	/* Free thresholded coefficients */
	free(preCompContrib -> compressedCoeffs);

	/* Free mRGBE encoded coefficients */
	free(preCompContrib -> mrgbeCoeffs);
	free(preCompContrib -> waveletFname);

	if (preCompContrib -> cache) {
	/* Free contribution cache */
	OOC_DeleteCache(preCompContrib -> cache);
	free(preCompContrib -> cache);
	}
	}

	/* Clean up precomputed contrib photon map */
	deletePhotons(preCompContribPmap);
	free(preCompContribPmap);
	}
	}



	void initPreComputedContrib (PreComputedContrib *preCompContrib)
	/* Initialise precomputed contribution container in photon map */
	{
	preCompContrib -> waveletFname = NULL;
	preCompContrib -> waveletFile = NULL;
	preCompContrib -> waveletMatrix =
	preCompContrib -> tWaveletMatrix = NULL;
	preCompContrib -> compressedCoeffs = NULL;
	setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0);
	setcolor(preCompContrib -> mrgbeRange.max, 0, 0, 0);
	setcolor(preCompContrib -> mrgbeRange.norm, 0, 0, 0);
	preCompContrib -> mrgbeCoeffs = NULL;
	preCompContrib -> scDim =
	preCompContrib -> nBins =
	preCompContrib -> coeffDim =
	preCompContrib -> nCoeffs =
	preCompContrib -> nNonZeroCoeffs =
	preCompContrib -> nCompressedCoeffs =
	preCompContrib -> contribSize = 0;
	preCompContrib -> cache = NULL;
	}



	void preComputeContrib (PhotonMap *pmap)
	/* Precompute contributions for a random subset of (finalGather *
	pmap -> numPhotons) photons, and return the per-modifier precomputed
	contribution photon maps in LUT preCompContribTab, discarding the
	original photons. */
	{
	unsigned long p, numPreComp;
	unsigned i, j, scDim, nBins,
	coeffDim, nCoeffs, nCompressedCoeffs,
	nNzDetailCoeffs;
	PhotonIdx pIdx;
	Photon photon;
	RAY ray;
	MODCONT *binnedContribs;
	COLR rgbe32 [WAVELET_PADD4_NUMAPPROX + 2];
	LUENT *preCompContribNode;
	PhotonMap *preCompContribPmap, nuPmap;
	PreComputedContrib *preCompContrib;
	LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode);
	WaveletMatrix2 waveletMatrix;
	#ifdef PMAP_CONTRIB_DBG
	RAY dbgRay;
	#endif

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

	/* Init new parent photon map and set output filename */
	initPhotonMap(&nuPmap, pmap -> type);
	nuPmap.fileName = pmap -> fileName;

	/* Set contrib/coeff mode in new parent photon map */
	nuPmap.contribMode = pmap -> contribMode;

	/* Allocate and init LUT containing per-modifier child photon maps */
	if (!(nuPmap.preCompContribTab = malloc(sizeof(LUTAB))))
	error(SYSTEM,
	"out of memory allocating LUT in preComputeContrib()"
	);
	memcpy(nuPmap.preCompContribTab, &lutInit, sizeof(LUTAB));

	/* Record start time, baby */
	repStartTime = time(NULL);
	#ifdef SIGCONT
	signal(SIGCONT, pmapPreCompReport);
	#endif
	repComplete = numPreComp = finalGather * pmap -> numPhotons;
	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 */
	pIdx = firstPhoton(pmap) + (unsigned long)(
	(p + 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 (i = 0; i < 3; i++)
	ray.ron [i] = photon.norm [i] / 127.0;


	/* Get contributions at photon position; retry with another
	photon if no contribs found */
	binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol);
	#ifdef PMAP_CONTRIB_DBG
	#if 0
	/* Get all contribs from same photon for debugging */
	/* Position must differ otherwise too many identical photon keys
	* will cause ooC octree to overflow! */
	VCOPY(dbgRay.rop, photon.pos);
	getPhoton(pmap, 0, &photon);
	memcpy(&dbgRay, &ray, sizeof(RAY));
	dbgRay.rsrc = photonSrcIdx(pmap, &photon);
	for (i = 0; i < 3; i++)
	dbgRay.ron [i] = photon.norm [i] / 127.0;
	binnedContribs = getPhotonContrib(pmap, &dbgRay, dbgRay.rcol);
	#endif
	#endif

	if (binnedContribs) {
	#if 0
	if (!(p & PMAP_CNTUPDATE))
	printf("Precomputing contribution photon %lu / %lu\n",
	p, numPreComp
	);
	#endif
	/* Found contributions at photon position, so generate
	* precomputed photon */
	if (!(preCompContribNode = lu_find(nuPmap.preCompContribTab,
	binnedContribs -> modname)
	)
	)
	error(SYSTEM, "out of memory allocating LUT entry in "
	"preComputeContrib()"
	);

	if (!preCompContribNode -> key) {
	/* New LUT node for precomputed contribs for this modifier */
	preCompContribNode -> key = (char*)binnedContribs -> modname;
	preCompContribNode -> data = (char*)(
	preCompContribPmap = malloc(sizeof(PhotonMap))
	);

	if (preCompContribPmap) {
	/* Init new child photon map and its contributions */
	initPhotonMap(preCompContribPmap,
	PMAP_TYPE_CONTRIB_CHILD
	);
	initPhotonHeap(preCompContribPmap);
	initContribHeap(preCompContribPmap);
	preCompContrib = (PreComputedContrib*)(
	preCompContribPmap -> preCompContrib =
	malloc(sizeof(PreComputedContrib))
	);
	initPreComputedContrib(preCompContrib);
	}
	if (!preCompContribPmap \|\| !preCompContrib)
	error(SYSTEM, "out of memory allocating new photon map "
	"in preComputeContrib()"
	);

	/* Set output filename from parent photon map */
	preCompContribPmap -> fileName = nuPmap.fileName;

	/* Set contrib/coeff mode from parent photon map */
	preCompContribPmap -> contribMode = nuPmap.contribMode;

	/* Set Shirley-Chiu square & wavelet matrix dimensions
	(number of bins resp. coefficients). */
	preCompContrib -> scDim = scDim = sqrt(
	preCompContrib -> nBins = nBins = binnedContribs -> nbins
	);

	if (scDim * scDim != nBins)
	/* Mkpmap shoulda took care of dis */
	error(INTERNAL, "contribution bin count not "
	"integer square in preComputeContrib()"
	);

	if (nBins > 1) {
	/* BINNING ENABLED; set up wavelet xform & compression.
	Get dimensions of wavelet coefficient matrix after
	boundary padding, and number of nonzero coefficients.
	The number of compressed (detail) coefficients is fixed
	and based on the number of NONZERO coefficients (minus
	approximations, see NOTE below) since zero coeffs are
	considered thresholded by default.

	!!! NOTE: THE APPROXIMATION COEFFICIENTS IN THE UPPER
	!!! LEFT OF THE MATRIX ARE _NEVER_ THRESHOLDED TO
	!!! MINIMISE COMPRESSION ARTEFACTS! THEY ARE ESSENTIAL
	!!! FOR RECONSTRUCTING THE ORIGINAL CONTRIBUTIONS. */
	preCompContrib -> coeffDim = coeffDim = padWaveletXform2(
	NULL, NULL, scDim, &preCompContrib -> nNonZeroCoeffs
	);
	preCompContrib -> nCoeffs = nCoeffs = coeffDim * coeffDim;
	nNzDetailCoeffs = WAVELET_PADD4_NUMDETAIL(
	preCompContrib -> nNonZeroCoeffs
	);
	nCompressedCoeffs = (1 - compressRatio) * nNzDetailCoeffs;

	if (nCompressedCoeffs < 1) {
	error(WARNING,
	"maximum contribution compression, clamping number "
	"of wavelet coefficients in preComputeContrib()"
	);
	nCompressedCoeffs = 1;
	}
	if (nCompressedCoeffs >= preCompContrib -> nCoeffs) {
	error(WARNING,
	"minimum contribution compression, clamping number "
	"of wavelet coefficients in preComputeContrib()"
	);
	nCompressedCoeffs = nNzDetailCoeffs;
	}
	preCompContrib -> nCompressedCoeffs = nCompressedCoeffs;

	/* Lazily allocate primary and transposed wavelet
	* coefficient matrix */
	preCompContrib -> waveletMatrix = waveletMatrix =
	allocWaveletMatrix2(coeffDim);
	preCompContrib -> tWaveletMatrix =
	allocWaveletMatrix2(coeffDim);
	if (!waveletMatrix \|\| !preCompContrib -> tWaveletMatrix)
	error(SYSTEM, "out of memory allocating wavelet "
	"coefficient matrix in preComputeContrib()"
	);

	/* Lazily allocate thresholded detail coefficients */
	preCompContrib -> compressedCoeffs = calloc(
	nNzDetailCoeffs, sizeof(PreComputedContribCoeff)
	);

	/* Lazily alloc mRGBE-encoded compressed wavelet coeffs */
	preCompContrib -> mrgbeCoeffs = calloc(
	nCompressedCoeffs, sizeof(mRGBE)
	);
	if (!preCompContrib -> compressedCoeffs \|\|
	!preCompContrib -> mrgbeCoeffs
	)
	error(SYSTEM, "out of memory allocating compressed "
	"contribution buffer in preComputeContrib()"
	);

	/* Set size of compressed contributions, in bytes */
	preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE(
	nCompressedCoeffs
	);
	}
	}
	else {
	/* LUT node already exists for this modifier; use its data */
	preCompContribPmap = (PhotonMap*)preCompContribNode -> data;
	preCompContrib = (PreComputedContrib*)(
	preCompContribPmap -> preCompContrib
	);
	scDim = preCompContrib -> scDim;
	nBins = preCompContrib -> nBins;
	nCoeffs = preCompContrib -> nCoeffs;
	nCompressedCoeffs = preCompContrib -> nCompressedCoeffs;
	waveletMatrix = preCompContrib -> waveletMatrix;
	}

	if (nBins > 1) {
	/* Binning enabled; copy binned contribs row-by-row to
	* wavelet xform input matrix, then compress & encode */
	for (i = 0; i < scDim; i++)
	memcpy(waveletMatrix [i],
	binnedContribs->cbin + PMAP_CONTRIB_XY2LIN(scDim, i, 0),
	scDim * WAVELET_COEFFSIZE
	);

	if (encodeContribs(preCompContrib, compressRatio))
	error(INTERNAL, "failed contribution "
	"compression/encoding in preComputeContrib()"
	);

	/* Encode wavelet approx coeffs in the upper left of the
	* wavelet matrix to 32-bit RGBE. These are not thresholded
	* to minimise compression artefacts. */
	for (i = 0; i < WAVELET_PADD4_APPROXDIM; i++)
	for (j = 0; j < WAVELET_PADD4_APPROXDIM; j++)
	setcolr(
	rgbe32 [PMAP_CONTRIB_XY2LIN(
	WAVELET_PADD4_APPROXDIM, i, j
	)],
	waveletMatrix [i][j][0], waveletMatrix [i][j][1],
	waveletMatrix [i][j][2]
	);

	/* Encode wavelet detail coeff range to 32-bit RGBE */
	setcolr(rgbe32 [WAVELET_PADD4_NUMAPPROX],
	preCompContrib -> mrgbeRange.max [0],
	preCompContrib -> mrgbeRange.max [1],
	preCompContrib -> mrgbeRange.max [2]
	);
	setcolr(rgbe32 [WAVELET_PADD4_NUMAPPROX + 1],
	preCompContrib -> mrgbeRange.min [0],
	preCompContrib -> mrgbeRange.min [1],
	preCompContrib -> mrgbeRange.min [2]
	);

	/* Dump 32-bit RGBE approx coeffs followed by mRGBE range to
	* contribution heap file. NOTE: mRGBE range minimum and
	* maximum are reversed in the file to facilitate handling
	* the special (but pointless) case of a single compressed
	* coeff; if so, only the mRGBE maximum is dumped, since the
	* minimum is implicitly zero and can be omitted to save
	* space */
	if (putbinary(rgbe32, sizeof(COLR),
	WAVELET_PADD4_NUMAPPROX + 1 + (nCompressedCoeffs > 1),
	preCompContribPmap -> contribHeap
	) != WAVELET_PADD4_NUMAPPROX + 1 + (nCompressedCoeffs > 1)
	)
	error(SYSTEM,
	"can't write wavelet approximation coefficients to "
	"contribution heap in preComputeContrib()"
	);

	/* Now dump mRGBE encoded, compressed detail coefficients */
	if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE),
	nCompressedCoeffs, preCompContribPmap -> contribHeap
	) != nCompressedCoeffs
	)
	error(SYSTEM,
	"can't write wavelet detail coefficients to "
	"contribution heap in preComputeContrib()"
	);
	}
	#if 1
	/* HACK: signal newPhoton() to set precomputed photon's
	contribution source from ray -> rsrc */
	/* XXX: DO WE STILL NEED THIS CRAP AFTER PRECOMP,
	SINCE CHILD PMAPS ARE PER-MODIFIER ANYWAY? */
	preCompContribPmap -> lastContribSrc.srcIdx = -2;
	#endif
	/* Append photon to new heap from ray and increment total
	count for all modifiers in parent photon map */
	newPhoton(preCompContribPmap, &ray);
	nuPmap.numPhotons++;
	}

	/* Update progress */
	repProgress++;

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

	/* Trash original pmap and its leaf file, and replace with new parent,
	which now acts as container for the per-modifier child pmaps */
	unlink(pmap -> store.leafFname);
	deletePhotons(pmap);
	memcpy(pmap, &nuPmap, sizeof(PhotonMap));

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

	/* Bail out if no photons could be precomputed */
	if (!pmap -> numPhotons)
	error(USER,
	"no contribution photons precomputed; try increasing -am"
	);

	/* Build per-modifier precomputed photon maps from their contribution
	heaps */
	lu_doall(pmap -> preCompContribTab, buildPreCompContribPmap, NULL);
	}



	#else
	/* -------------------------------------------------------------------
	U N I T T E S T S
	------------------------------------------------------------------- */

	#include <stdio.h>
	#include <random.h>
	#include "func.h"

	int main (int argc, char *argv [])
	{
	unsigned i, scdim, nsamp, numTheta = 0, numPhi = 0;
	double t, p;
	RAY ray;
	int bin;

	if (argc < 3) {
	fprintf(stderr, "%s <scdim> <nsamp> [<var>=<value>; ..]\n",
	argv [0]
	);
	fprintf(stderr, "Default variable defs: %s\n",
	PMAP_CONTRIB_SCDEFAULTS
	);
	fputs("\nMissing Shirley-Chiu dimension scDim>1, "
	"number of samples nsamp\n",
	stderr
	);
	return -1;
	}

	if (!(scdim = atoi(argv [1]))) {
	fputs("Invalid Shirley-Chiu dimension\n", stderr);
	return -1;
	}

	if (!(nsamp = atoi(argv [2]))) {
	fputs("Invalid num samples\n", stderr);
	return -1;
	}
	else {
	numTheta = (int)(sqrt(nsamp) / 2);
	numPhi = 4* numTheta;
	}

	/* (Doan' need to) Init cal func routines for binning */
	#if 0
	initfunc();
	setcontext(RCCONTEXT);
	#endif

	/* Compile default orientation variables for contribution binning */
	scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0);
	/* Compile custom orientation variabls from command line */
	for (i = 3; i < argc; i++)
	scompile(argv [i], NULL, 0);

	for (i = 0; i < nsamp; i++) {
	#if 0
	/* Random */
	t = 0.5 * PI * drand48();
	p = 2 * PI * drand48();
	#else
	/* Stratified */
	t = 0.5 * PI * ((i % numTheta) + drand48()) / numTheta;
	p = 2.0 * PI * ((i / numTheta) + drand48()) / numPhi;
	#endif
	ray.rdir [0] = sin(t) * cos(p);
	ray.rdir [1] = sin(t) * sin(p);
	ray.rdir [2] = cos(t);
	bin = ray2bin(&ray, scdim);

	printf("%.3f\t%.3f\t%.3f\t-->\t%d\n",
	ray.rdir [0], ray.rdir [1], ray.rdir [2], bin
	);
	}

	return 0;
	}
	#endif

	#endif /* PMAP_CONTRIB */

pmapcontrib.cNo OneTemporaryActions

File Metadata

pmapcontrib.cView Options

Event Timeline

pmapcontrib.c
No OneTemporary
Actions

pmapcontrib.c
View Options