Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmapcontrib.c
45f043963340pmapmerge

Last Change
Raw File

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 $
*/
/* TODO
-- Handle 1 bin (i.e. skip wavelet compression)
-- Save per-mod photon maps
*/
#include "pmapcontrib.h"
#include "pmapdiag.h"
#include "pmaprand.h"
#include "pmapmat.h"
#include "pmapsrc.h"
#include "pmapio.h"
#include "pmutil.h"
#include "otspecial.h"
#include "otypes.h"
#include "lookup.h"
extern void SDdisk2square(double sq[2], double diskx, double disky);
static int logDim (unsigned size)
/* Return log2(sqrt(size)) = l, where size = (2^l)(2^l).
Is there a faster way (e.g. binary search)? */
{
unsigned i, sz, dim;
if (!size)
return 0;
for (i = 0, dim = 0, sz = size; sz >>= 2; dim++);
return dim;
}
static int xy2bin (unsigned l, int x, int y)
/* Serialise 2D coordinates in range (2^l) x (2^l) to 1D bin.
Returns -1 if coordinates invalid */
{
return x < 0 || y < 0
? -1
: (x << l) + y;
}
static void bin2xy (unsigned l, int bin, int *x, int *y)
/* Deserialise 1D bin to 2D coordinates in range (2^l) x (2^l).
Returns -1 if bin invalid */
{
/* Obviously this is faster than integer division/modulo */
*x = bin < 0 ? -1 : bin >> l;
*y = bin < 0 ? -1 : bin & PMAP_CONTRIB_SCDIM(l) - 1;
}
static int ray2bin (const RAY *ray, unsigned nbins)
/* Map ray dir (pointing away from origin) to its 1D Shirley-Chiu bin,
where nbins = (2^l) x (2^l) for l > 1.
Returns -1 if mapped bin is invalid (e.g. behind plane) */
{
const unsigned l = logDim(nbins), scDim = PMAP_CONTRIB_SCDIM(l);
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 - dz*dz) / 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 xy2bin(l, scBin [0], scBin [1]);
}
else return -1;
}
#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;
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;
worldfunc(RCCONTEXT, &srcRay);
set_eparams((char*)srcCont -> params);
if ((bin = ray2bin(&srcRay, 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
);
}
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(
pmap -> contribTab, 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;
}
static int coeffCompare (const void *c1, const void *c2)
/* Comparison function to REVERSE sort thresholded coefficients */
{
const ThreshWaveletCoeff *tcoeff1 = c1,
*tcoeff2 = c2;
/* TODO: USE MAX HERE INSTEAD OF AVERAGE? */
const WAVELET_COEFF v1 = colorAvg(tcoeff1 -> coeff),
v2 = colorAvg(tcoeff2 -> coeff);
if (v1 < v2)
return 1;
else if (v1 > v2)
return -1;
else
return 0;
}
static int thresholdContribs (MODCONT *binnedContribs,
PreComputedContrib *preCompContrib
)
/* Threshold binned detail coefficients in binnedContribs -> cbin [1]
by keeping the (preCompContrib -> nCompressedBins) largest and
returning these in preContrib -> threshCoeffs along with their
original bin order.
NOTE: binnedContribs -> cbin [0] is the average wavelet coefficient
and excluded from thresholding.
Returns 0 on success.
*/
{
unsigned b;
ThreshWaveletCoeff *threshCoeffs;
if (!preCompContrib)
return -1;
/* Skip 1st coeff; it's the average and therefore not thresholded */
for (b = 0, threshCoeffs = preCompContrib -> threshCoeffs;
b < binnedContribs -> nbins - 1;
b++
) {
/* Set up pointers to coeffs (sorted faster than 3 floats/doubles)
and remember original bin position prior to sorting.
NOTE: The original bin position *excludes* the 1st coeff,
i.e. the average! */
threshCoeffs [b].coeff = (WAVELET_COEFF*)(
&(binnedContribs->cbin [b+1])
);
threshCoeffs [b].bin = b;
}
/* REVERSE sort coeffs by magnitude; the non-thresholded coeffs
will then start at threshCoeffs [0] */
qsort(threshCoeffs, binnedContribs -> nbins - 1,
sizeof(ThreshWaveletCoeff), coeffCompare
);
return 0;
}
static int encodeContribs (MODCONT *binnedContribs,
PreComputedContrib *preCompContrib, float compressRatio
)
/* Apply wavelet transform to binnedContribs -> cbin and compress
according to compressRatio, storing thresholded and mRGBE-encoded
coefficients in preCompContrib -> mrgbeCoeffs.
Note that the average coefficient is not encoded, and returned
in binnedContribs -> cbin [0].
Returns 0 on success. */
{
unsigned b, i;
ThreshWaveletCoeff *threshCoeffs;
mRGBERange *mrgbeRange;
mRGBE *mrgbeCoeff;
if (!binnedContribs || !preCompContrib)
return -1;
/* Do 2D wavelet transform on preCompContrib -> waveletMatrix
(which really maps to binnedContribs -> cbin) */
if (waveletXform2(preCompContrib -> waveletMatrix,
preCompContrib -> tWaveletMatrix, preCompContrib -> l
) < 0
)
error(INTERNAL, "failed wavelet transform of contributions");
/* Compress wavelet detail coeffs by thresholding */
if (thresholdContribs(binnedContribs, preCompContrib) < 0)
error(INTERNAL, "failed wavelet compression of contributions");
threshCoeffs = preCompContrib -> threshCoeffs;
/* 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 (b = 0; b < preCompContrib -> nCompressedBins; b++) {
for (i = 0; i < 3; i++) {
if (threshCoeffs [b].coeff [i] < mrgbeRange -> min [i])
mrgbeRange -> min [i] = threshCoeffs [b].coeff [i];
if (threshCoeffs [b].coeff [i] > mrgbeRange -> max [i])
mrgbeRange -> max [i] = threshCoeffs [b].coeff [i];
}
}
/* Init mRGBE coefficient normalisation from range */
mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max);
mrgbeCoeff = preCompContrib -> mrgbeCoeffs;
/* Encode wavelet detail coefficients to mRGBE */
for (b = 0; b < preCompContrib -> nCompressedBins; b++)
mrgbeCoeff [b] = mRGBEencode(threshCoeffs [b].coeff,
mrgbeRange, threshCoeffs [b].bin
);
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(
pmap -> contribTab, 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 photon distance */
scalecolor(flux, 1 - sqn -> dist2 / r2);
#endif
addcolor(irrad, flux);
addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux);
}
return contrib;
}
static void freePreCompContribNode (void *p)
/* Free per-modifier precomputed contributions LUT entry */
{
PreComputedContrib *preCompContrib = p;
/* Free primary and transposed wavelet matrices */
free(preCompContrib -> waveletMatrix);
freeWaveletMatrix(preCompContrib -> tWaveletMatrix,
preCompContrib -> l
);
/* Free thresholded coefficients */
free(preCompContrib -> threshCoeffs);
/* Free mRGBE encoded coefficients */
free(preCompContrib -> mrgbeCoeffs);
/* Clean up precomputed contrib photon map */
if (preCompContrib -> pmap) {
deletePhotons(preCompContrib -> pmap);
free(preCompContrib -> pmap);
}
}
int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p)
/* Build per-modifier precomputed photon map from compressed
contributions in preCompContrib -> mrgbeCoeffs and save to a common
subdirectory along with the compressed wavelet coefficients.
Returns 0 on success. */
{
PreComputedContrib *preCompContrib = (PreComputedContrib*)(
preCompContribNode -> data
);
PhotonMap *preCompContribPmap = preCompContrib -> pmap;
/* Flush photon/contrib heaps */
flushPhotonHeap(preCompContribPmap);
fflush(preCompContribPmap -> contribHeap);
if (verbose) {
sprintf(errmsg, "Building precomputed contribution photon map "
"for modifier %s\n", preCompContribNode -> key
);
eputs(errmsg);
#if NIX
fflush(stderr);
#endif
}
/* Rebuild underlying data structure, destroying heap */
buildPhotonMap(preCompContribPmap, NULL, NULL, 1);
/* Free primary and transposed wavelet matrices; no longer needed */
free(preCompContrib -> waveletMatrix);
freeWaveletMatrix(preCompContrib -> tWaveletMatrix, preCompContrib -> l);
/* Free thresholded coefficients; no longer needed */
free(preCompContrib -> threshCoeffs);
preCompContrib -> waveletMatrix =
preCompContrib -> tWaveletMatrix = NULL;
preCompContrib -> threshCoeffs = NULL;
return 0;
}
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. */
/* !!! NOTE: PRECOMPUTATION CURRENTLY WITHOUT OOC CACHE !!! */
{
unsigned long i, numPreComp;
unsigned b, j;
long nCompressedBins;
PhotonIdx pIdx;
Photon photon;
RAY ray;
MODCONT *binnedContribs;
COLR coeffNorm;
LUENT *preCompContribNode;
PreComputedContrib *preCompContribData;
PhotonMap *preCompContribPmap, nuPmap;
LUTAB lutInit = LU_SINIT(NULL,
freePreCompContribNode
);
if (verbose) {
sprintf(errmsg,
"\nPrecomputing contributions for %ld photons\n", numPreComp
);
eputs(errmsg);
#if NIX
fflush(stderr);
#endif
}
/* Init new photon map and set output filename */
initPhotonMap(&nuPmap, pmap -> type);
nuPmap.fileName = pmap -> fileName;
/* Allocate and init LUT containing per-modifier precomputed
* contribution 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 (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 = getPhotonContrib(pmap, &ray, ray.rcol);
} while (!binnedContribs);
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;
preCompContribData = (PreComputedContrib*)(
preCompContribNode -> data =
malloc(sizeof(PreComputedContrib))
);
if (preCompContribData)
preCompContribPmap = preCompContribData -> pmap =
malloc(sizeof(PhotonMap));
if (!preCompContribData || !preCompContribPmap)
error(SYSTEM, "out of memory allocating photon map in "
"preComputeContrib()"
);
/* Init new precomputed photon map and contribution heap */
initPhotonMap(preCompContribPmap, nuPmap.type);
initPhotonHeap(preCompContribPmap);
initContribHeap(preCompContribPmap);
/* Get Shirley-Chiu square & wavelet matrix resolution, bins
per dimension, and fixed number of compressed coeffs/bins
(note we subtract one as the average coefficient is not
thresholded) */
preCompContribData -> l = logDim(binnedContribs -> nbins);
preCompContribData -> scDim =
PMAP_CONTRIB_SCDIM(preCompContribData -> l);
if (binnedContribs -> nbins > 1) {
/* Binning enabled; set up wavelet xform & compression */
nCompressedBins =
(binnedContribs -> nbins - 1) * (1 - compressRatio);
if (nCompressedBins > 0)
preCompContribData -> nCompressedBins = nCompressedBins;
else error(USER,
"invalid compression ratio in preComputeContrib()"
);
/* Lazily "allocate" primary wavelet matrix pointing to array
binnedContribs -> cbin; this imposes the 2D matrix structure
expected by the wavelet xform routines onto the 1D array. */
if (!(preCompContribData -> waveletMatrix = calloc(
preCompContribData -> scDim, sizeof(WaveletCoeff3*))
)
)
error(SYSTEM, "out of memory allocating primary wavelet "
"matrix in preComputeContrib()"
);
for (b = 0; b < preCompContribData -> scDim; b++)
/* Point to each row in existing 1D contrib array */
preCompContribData -> waveletMatrix [b] =
&binnedContribs -> cbin [b * preCompContribData -> scDim];
/* Lazily allocate transposed wavelet matrix */
preCompContribData -> tWaveletMatrix =
allocWaveletMatrix(preCompContribData -> l);
if (!preCompContribData -> tWaveletMatrix)
error(SYSTEM, "out of memory allocating transposed wavelet "
"matrix in preComputeContrib()"
);
/* Lazily allocate thresholded detail coefficients (minus
the average coeff) */
preCompContribData -> threshCoeffs = calloc(
binnedContribs -> nbins - 1, sizeof(ThreshWaveletCoeff)
);
/* Lazily allocate mRGBE-encoded, compressed wavelet coeffs */
preCompContribData -> mrgbeCoeffs = calloc(
preCompContribData -> nCompressedBins, sizeof(mRGBE)
);
if (!preCompContribData -> threshCoeffs ||
!preCompContribData -> mrgbeCoeffs
)
error(SYSTEM, "out of memory allocating compressed "
"coefficients in preComputeContrib()"
);
}
else preCompContribData =
(PreComputedContrib*)preCompContribNode -> data;
/* Compress & encode binned contribs */
if (encodeContribs(binnedContribs, preCompContribData,
compressRatio
) < 0
)
error(INTERNAL, "failed contribution compression/encoding "
"in preComputeContrib()"
);
/* Dump encoded bins to contrib heap file, prepended by
normalisation in 32-bit RGBE */
setcolr(coeffNorm,
preCompContribData -> mrgbeRange.norm [0],
preCompContribData -> mrgbeRange.norm [1],
preCompContribData -> mrgbeRange.norm [2]
);
putbinary(coeffNorm, 1, sizeof(coeffNorm),
preCompContribPmap -> contribHeap
);
if (putbinary(preCompContribData -> mrgbeCoeffs, sizeof(mRGBE),
preCompContribData -> nCompressedBins,
preCompContribPmap -> contribHeap
) != preCompContribData -> nCompressedBins
)
error(SYSTEM, "failed writing to coefficients to "
"contribution heap in preComputeContrib()"
);
}
else {
/* No binning; skip wavelet xform & compression */
preCompContribData -> nCompressedBins = 0;
preCompContribData -> waveletMatrix =
preCompContribData -> tWaveletMatrix = NULL;
preCompContribData -> threshCoeffs = NULL;
preCompContribData -> mrgbeCoeffs = NULL;
}
/* Set photon flux to coarsest average wavelet coefficient
NOTE: binnedContrib -> cbin [0] ==
preCompContribData -> waveletMatrix [0][0] ==
preCompContribData -> tWaveletMatrix [0][0].
If binning is disabled, this is the total contribution from all
directions */
copycolor(ray.rcol, binnedContribs -> cbin [0]);
/* HACK: signal newPhoton() to set precomputed photon's
contribution source from ray -> rsrc */
preCompContribPmap -> lastContribSrc.srcIdx = -2;
/* 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 replace with new one, which now acts as a
* container for the constituent per-modifier precomputed pmaps */
deletePhotons(pmap);
memcpy(pmap, &nuPmap, sizeof(PhotonMap));
#ifdef SIGCONT
signal(SIGCONT, SIG_DFL);
#endif
/* 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, l, nbins, nsamp, numTheta = 0, numPhi = 0;
double t, p;
RAY ray;
int bin;
if (argc < 3) {
fprintf(stderr, "%s <l> <nsamp> [<var>=<value>; ..]\n", argv [0]);
fprintf(stderr, "Default variable defs: %s\n",
PMAP_CONTRIB_SCDEFAULTS
);
fputs("\nMissing resolution l>1, number of samples nsamp\n",
stderr
);
return -1;
}
if (!(l = atoi(argv [1]))) {
fputs("Invalid resolution\n", stderr);
return -1;
}
else nbins = PMAP_CONTRIB_SCBINS(l);
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, nbins);
printf("%.3f\t%.3f\t%.3f\t-->\t%d\n",
ray.rdir [0], ray.rdir [1], ray.rdir [2], bin
);
}
return 0;
}
#endif
c4science · Help