Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmcontrib4.c
9f35cb983d38pmapmerge

Last Change
Raw File

pmcontrib4.c
View Options

#ifndef lint
static const char RCSid[] = "$Id: pmcontrib2.c,v 2.5 2018/11/08 00:54:07 greg Exp $";
#endif
/*
=========================================================================
Photon map routines for precomputed light source contributions.
These routines interface to rcontrib.
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$
*/
#include "pmcontribcache.h"
#ifdef PMAP_CONTRIB
#include "pmapmat.h"
#include "pmapsrc.h"
#include "pmaprand.h"
#include "pmapio.h"
#include "pmapdiag.h"
#include "otypes.h"
#include "otspecial.h"
#include "oocnn.h"
#include "random.h"
void initPmapContribTab (LUTAB *contribTab, int *contribMode)
/* Set contribution table (interface to rcmain.c) */
{
pmapContribTab = contribTab;
pmapContribMode = contribMode;
}
/* ------------------ CONTRIB PMAP LOADING STUFF --------------------- */
static int checkContribModifier (const LUENT *modContNode, void *p)
/* Check for valid modifier in LUT entry */
{
MODCONT *modCont = (MODCONT*)modContNode -> data;
char *modName = (char*)modCont -> modname;
OBJECT modObj = modifier(modName);
const PhotonMap *pmap = (PhotonMap*)p;
if (modObj == OVOID) {
sprintf(errmsg, "invalid modifier %s", modName);
error(USER, errmsg);
}
if (!islight(objptr(modObj) -> otype)) {
sprintf(errmsg, "%s is not a light source modifier", modName);
error(USER, errmsg);
}
/* Warn if number of bins exceeds lookup bandwidth, as this is
* guaranteed to result in empty bins! */
if (pmap -> maxGather < modCont -> nbins) {
sprintf(errmsg, "contribution photon lookup bandwidth too low "
"for modifier %s, some bins will be zero", modName
);
error(WARNING, errmsg);
}
return 0;
}
void initPmapContrib (PhotonMap *pmap)
/* Set up photon map contributions and get binning parameters */
{
unsigned t;
/* Avoid incompatible photon maps */
for (t = 0; t < NUM_PMAP_TYPES; t++)
if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) {
sprintf(errmsg, "%s photon map does not support contributions",
pmapName [t]
);
error(USER, errmsg);
}
if (!pmapContribTab)
error(INTERNAL,
"contribution LUT not initialised in initPmapContrib()"
);
/* Check for valid contribution modifiers */
lu_doall(pmapContribTab, checkContribModifier, pmap);
/* Set callback to collect contributions */
if (pmap) {
pmap -> lookup = getPreCompPhotonContrib;
/* XXX: Need tighter tolerance for contribution photon lookups? */
pmap -> gatherTolerance = 1.0;
}
}
static int loadPreCompContrib (const LUENT *modContNode, void *p)
/* Load child contribution photon map for current contrib modifier entry */
{
MODCONT *modCont = (MODCONT*)modContNode -> data;
char *modName = (char*)modCont -> modname;
PhotonMap *preCompContribPmap;
const PhotonMap *parentPmap = (const PhotonMap*)p;
LUENT *preCompContribNode;
PreComputedContrib *preCompContrib;
char dirName [1024];
/* Allocate new LUT entry for child precomputed contribution pmap */
preCompContribNode = lu_find(parentPmap -> preCompContribTab, modName);
if (!preCompContribNode)
error(SYSTEM,
"out of memory allocating LUT entry in loadPreCompContrib()"
);
preCompContribNode -> key = modName;
/* Allocate child precomputed contribution photon map */
preCompContribNode -> data = (char*)(
preCompContribPmap = malloc(sizeof(PhotonMap))
);
if (!preCompContribPmap)
error(SYSTEM, "out of memory allocating precomputed contribution "
"photon map in loadPreCompContrib()"
);
/* Set subdirectory from parent photon map's filename */
snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR,
parentPmap -> fileName
);
/* Allocate & set child pmap's filename from subdir and modifier */
preCompContribPmap -> fileName = malloc(strlen(parentPmap -> fileName) +
strlen(modName) + strlen(PMAP_CONTRIB_FILE)
);
if (!preCompContribPmap -> fileName)
error(SYSTEM, "out of memory allocating filename in "
"loadPreCompContrib()"
);
snprintf(preCompContribPmap -> fileName, sizeof(dirName) + 4,
PMAP_CONTRIB_FILE, dirName, modName
);
/* Load child precomp contrib photon map for current modifier;
loadPhotonMap() will not recurse and call loadContribPhotonMap()
again because preCompContribPmap -> preCompContribTab == NULL */
loadPhotonMap(preCompContribPmap, preCompContribPmap -> fileName);
/* Pass parent pmap's lookup parameters on to child */
preCompContribPmap -> minGather = parentPmap -> minGather;
preCompContribPmap -> maxGather = parentPmap -> maxGather;
preCompContribPmap -> gatherTolerance = parentPmap -> gatherTolerance;
/* Override contrib/coefficient mode if it doesn't match precomputed */
if (*pmapContribMode != preCompContribPmap -> contribMode) {
sprintf(errmsg,
"photon map contains precomputed %s, overriding rcontrib setting",
preCompContribPmap -> contribMode ? "contributions" : "coefficients"
);
error(WARNING, errmsg);
*pmapContribMode = preCompContribPmap -> contribMode;
}
/* NOTE: preCompContribPmap -> preCompContrib is initialised by
* loadPhotonMap() */
preCompContrib = (PreComputedContrib*)(
preCompContribPmap -> preCompContrib
);
/* Set number of bins and wavelet coefficients */
preCompContrib -> nBins = preCompContrib -> scDim *
preCompContrib -> scDim;
preCompContrib -> nCoeffs = preCompContrib -> coeffDim *
preCompContrib -> coeffDim;
/* Set wavelet coefficient filename */
preCompContrib -> waveletFname = malloc(strlen(dirName) +
strlen(modName) + strlen(PMAP_CONTRIB_WAVELETFILE)
);
if (!preCompContribPmap -> fileName)
error(SYSTEM, "out of memory allocating filename in "
"loadPreCompContrib()"
);
snprintf(preCompContrib -> waveletFname, sizeof(dirName) + 5,
PMAP_CONTRIB_WAVELETFILE, dirName, modName
);
return 0;
}
void loadContribPhotonMap (PhotonMap *pmap, const char *fname)
/* Load contribution pmap and its per-modifier precomputed children */
{
LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode);
/* Allocate & init LUT containing per-modifier child contrib pmaps */
if (!(pmap -> preCompContribTab = malloc(sizeof(lutInit))))
error(SYSTEM, "out of memory allocating precomputed contribution "
"LUT in loadContribPhotonMap()"
);
memcpy(pmap -> preCompContribTab, &lutInit, sizeof(lutInit));
/* NOTE: filename already set in loadPhotonMap()
pmap -> fileName = (char*)fname; */
/* Init contrib photon map for light source contributions */
initPmapContrib(pmap);
/* Load child contribution photon maps for each modifier in contribTab */
lu_doall(pmapContribTab, loadPreCompContrib, pmap);
}
static int decodeContribs (PreComputedContrib *preCompContrib)
/* Decode and decompress mRGBE-encoded wavelet coefficients in
preCompContrib -> mrgbeCoeffs, apply inverse wavelet transform and
return decoded contributions in the wavelet coefficient matrix
preCompContrib -> waveletMatrix.
NOTE: THE WAVELET COEFFICIENT MATRIX IS ASSUMED TO BE ZEROED, with
the average wavelet coefficient in the upper left of the matrix
(i.e. preCompContrib -> waveletMatrix [0][0]).
Returns 0 on success, else -1. */
{
unsigned c, i, j, coeffIdx, scDim,
coeffDim, nCoeffs, nCompressedCoeffs;
WaveletMatrix2 waveletMatrix;
mRGBERange *mrgbeRange;
mRGBE *mrgbeCoeffs;
DCOLOR fCoeff;
if (!preCompContrib || !preCompContrib -> mrgbeCoeffs ||
!(waveletMatrix = preCompContrib -> waveletMatrix) ||
!preCompContrib -> tWaveletMatrix
)
/* Should be initialised by caller! */
return -1;
scDim = preCompContrib -> scDim;
coeffDim = preCompContrib -> coeffDim;
nCoeffs = preCompContrib -> nCoeffs;
nCompressedCoeffs = preCompContrib -> nCompressedCoeffs;
mrgbeCoeffs = preCompContrib -> mrgbeCoeffs;
mrgbeRange = &preCompContrib -> mrgbeRange;
/* Init mRGBE coefficient normalisation from range */
if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max))
return -1;
/* Decode mRGBE coefficients and populate wavelet coefficient matrix,
omitting average at waveletMatrix [0][0].
NOTE: The matrix should be initialised to zero so */
for (c = 0; c < nCompressedCoeffs; c++) {
coeffIdx = mRGBEdecode(mrgbeCoeffs [c], mrgbeRange, fCoeff);
#ifdef PMAP_CONTRIB_DBG
/* Check for invalid decoded coefficients */
if (coeffIdx < 1 || coeffIdx >= nCoeffs)
error(CONSISTENCY,
"wavelet coefficient index out of range in decodeContribs()"
);
#endif
i = PMAP_CONTRIB_LIN2X(coeffDim, coeffIdx);
j = PMAP_CONTRIB_LIN2Y(coeffDim, coeffIdx);
copycolor(waveletMatrix [i] [j], fCoeff);
}
/* Do inverse 2D wavelet transform on preCompContrib -> waveletMatrix
(which actually maps to binnedContribs) */
if (!padWaveletInvXform2(waveletMatrix,
preCompContrib -> tWaveletMatrix, coeffDim, scDim
)
)
error(INTERNAL,
"failed inverse wavelet transform in decodeContribs()"
);
#ifdef PMAP_CONTRIB_DBG
for (i = 0; i < scDim; i++) {
for (j = 0; j < scDim; j++) {
/* Dump decoded contribs for debugging */
printf("% 7.3g\t", colorAvg(waveletMatrix [i][j]));
/* Warn if coefficient is negative; this _can_ happen due to loss
of precision by the mRGBE encoding */
if (min(min(waveletMatrix [i][j][0], waveletMatrix [i][j][1]),
waveletMatrix [i][j][2]
) < 0
)
error(WARNING, "negative contributions in getPreCompContrib()");
}
putchar('\n');
}
putchar('\n');
#endif
return 0;
}
static void getPreCompContribByPhoton (const Photon *photon,
OOC_DataIdx photonIdx, PreComputedContrib *preCompContrib,
DCOLOR *binnedContribs
)
/* Fetch and decode mRGBE-encoded wavelet coefficients for given photon
from wavelet file at offset photonIdx, perform inverse wavelet xform,
and return the reconstructed binned contributions in binnedContribs
(which is assumed to be variable). Returns 0 on success, else -1. */
/* XXX: HERE binnedContribs POINTS TO CACHE PAGE DATA! */
{
unsigned i, j, scDim, coeffDim, nCompressedCoeffs;
COLR mrgbeRange32 [2];
COLOR fCoeff;
if (!binnedContribs || !preCompContrib)
/* Shouldn't happen */
error(INTERNAL,
"contributions not initialised in getPreCompContrib()"
);
if (preCompContrib -> nBins <= 1)
/* No binning, so nothing to decode */
return;
scDim = preCompContrib -> scDim;
coeffDim = preCompContrib -> coeffDim;
nCompressedCoeffs = preCompContrib -> nCompressedCoeffs;
/* Lazily initialise preCompContrib */
if (!preCompContrib -> waveletFile) {
/* Lazily open wavelet coefficient file */
preCompContrib -> waveletFile = fopen(
preCompContrib -> waveletFname, "rb"
);
if (!preCompContrib -> waveletFile) {
sprintf(errmsg, "failed to open wavelet coefficient file %s "
"in getPreCompContrib()", preCompContrib -> waveletFname
);
error(SYSTEM, errmsg);
}
/* Set record size of encoded coeffs in wavelet file */
preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE(
nCompressedCoeffs
);
/* Lazily allocate buffer for mRGBE wavelet coeffs */
preCompContrib -> mrgbeCoeffs = calloc(
nCompressedCoeffs, sizeof(mRGBE)
);
if (!preCompContrib -> mrgbeCoeffs)
error(SYSTEM, "out of memory allocating decoded wavelet "
"coefficients in getPreCompContrib()"
);
/* Lazily allocate primary and transposed wavelet matrices */
preCompContrib -> waveletMatrix = allocWaveletMatrix2(coeffDim);
preCompContrib -> tWaveletMatrix = allocWaveletMatrix2(coeffDim);
if (!preCompContrib -> waveletMatrix ||
!preCompContrib -> tWaveletMatrix
)
error(SYSTEM, "out of memory allocating wavelet coefficient "
"matrix in getPreCompContrib()"
);
}
/* Seek to photon index in wavelet file and read mRGBE normalisation
and encoded bins */
if (fseek(preCompContrib -> waveletFile,
photonIdx * preCompContrib -> contribSize, SEEK_SET
) < 0
)
error(SYSTEM, "failed seek in wavelet coefficient file "
"in getPreCompContrib()"
);
/* Read only mRGBE maximum if 1 compressed bin (maximum compression),
* since minimum is implicitly zero and therefore omitted */
getbinary(mrgbeRange32, sizeof(COLR), 1 + (nCompressedCoeffs > 1),
preCompContrib -> waveletFile
);
if (getbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE),
nCompressedCoeffs, preCompContrib -> waveletFile
) != nCompressedCoeffs
)
error(SYSTEM, "failed reading from wavelet coefficient file "
"in getPreCompContrib()"
);
/* Get mRGBE range (NOTE min/max are reversed in the coeff file) */
colr_color(fCoeff, mrgbeRange32 [0]);
copycolor(preCompContrib -> mrgbeRange.max, fCoeff);
if (nCompressedCoeffs > 1) {
colr_color(fCoeff, mrgbeRange32 [1]);
copycolor(preCompContrib -> mrgbeRange.min, fCoeff);
}
else {
/* Single compressed bin; mRGBE minimum is implicitly zero */
setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0);
}
/* Zero wavelet coefficient matrix and set average coefficient from
photon flux; this coefficient is expected by the inverse transform
routine in the upper left of the matrix (i.e. preCompContrib ->
waveletMatrix [0] [0]) */
zeroCoeffs2(preCompContrib -> waveletMatrix, coeffDim);
getPhotonFlux((Photon*)photon, fCoeff);
copycolor(preCompContrib -> waveletMatrix [0] [0], fCoeff);
/* All set, now decode mRGBE coeffs and invert wavelet transform */
if (decodeContribs(preCompContrib))
error(INTERNAL, "failed contribution decoding/decompression "
"in getPreCompContrib()"
);
/* Copy decoded contributions from wavelet coefficient matrix rows */
for (i = 0; i < scDim; i++)
memcpy(binnedContribs + PMAP_CONTRIB_XY2LIN(scDim, i, 0),
preCompContrib -> waveletMatrix [i], scDim * WAVELET_COEFFSIZE
);
}
typedef struct {
const RAY *ray;
RREAL rayCoeff [3];
COLORV *totalContrib;
} PreCompContribRCData;
static int getPreCompContribByMod (const LUENT *preCompContribTabNode,
void *p
)
/* Fetch and decode precomputed contributions from closest photon in pmap
* referenced by preCompContribTabNode, and accumulate in pmapContribTab
* for the corresponding modifier */
{
PhotonMap *preCompContribPmap = (PhotonMap*)(
preCompContribTabNode -> data
);
PreComputedContrib *preCompContrib;
const char *modName = preCompContribTabNode -> key;
MODCONT *modCont;
PreCompContribRCData *rcData = (PreCompContribRCData*)p;
LUENT *contribTabNode;
Photon photon;
OOC_DataIdx photonIdx;
COLOR photonContrib;
unsigned b, k;
if (!preCompContribPmap || !rcData ||
!(preCompContrib = (PreComputedContrib*)(
preCompContribPmap -> preCompContrib
))
) {
sprintf(errmsg, "missing precomputed contributions for modifier %s",
modName
);
error(INTERNAL, errmsg);
}
/* Find rcontrib's LUT entry for this modifier */
contribTabNode = lu_find(pmapContribTab, modName);
if (!contribTabNode || !contribTabNode -> key ||
!(modCont = (MODCONT*)contribTabNode -> data)
) {
sprintf(errmsg, "missing contribution LUT entry for modifier %s",
modName
);
error(INTERNAL, errmsg);
}
/* Reallocate rcontrib bins if they don't match precomputed */
if (modCont -> nbins != preCompContrib -> nBins) {
contribTabNode -> data = realloc(modCont,
sizeof(MODCONT) + sizeof(DCOLOR) * (preCompContrib -> nBins - 1)
);
if (!contribTabNode -> data) {
sprintf(errmsg, "out of memory reallocating contribution bins "
"for modifier %s", modCont -> modname
);
error(SYSTEM, errmsg);
}
else {
sprintf(errmsg, "bin count mismatch for modifier %s, resizing",
modCont -> modname
);
error(WARNING, errmsg);
}
modCont = (MODCONT*)contribTabNode -> data;
modCont -> nbins = preCompContrib -> nBins;
memset(modCont -> cbin, 0, sizeof(DCOLOR) * modCont -> nbins);
}
/* Fetch closest photon and its contributions */
if (find1Photon(preCompContribPmap, rcData -> ray, &photon,
&photonIdx
)
) {
if (preCompContrib -> nBins > 1) {
/* BINNING ENABLED; lazily init contrib cache if necessary */
if (initContribCache(preCompContrib)) {
/* CACHE ENABLED; fetch cache page for found photon */
DCOLOR *cacheBins;
if (!getContribCache(preCompContrib, photonIdx, &cacheBins)) {
/* PAGE NOT CACHED; decode contribs into new cache page */
getPreCompContribByPhoton(&photon, photonIdx, preCompContrib,
cacheBins
);
}
else; /* PAGE CACHED; (re)use decoded contribs from cache! */
for (b = 0; b < preCompContrib -> nBins; b++) {
/* Accumulate binned contributions into rcontrib's LUT entry
for this modifier, weighted by ray coefficient. Also sum
total contributions. */
/* NOTE: Using multcolor() on the decoded contribs would
modify them in the cache, so use a temp variable here. */
for (k = 0; k < 3; k++)
photonContrib [k] = cacheBins [b] [k] *
rcData -> rayCoeff [k];
addcolor(modCont -> cbin [b], photonContrib);
addcolor(rcData -> totalContrib, photonContrib);
}
}
else {
/* CACHE DISABLED; decode contribs into temp array on stack */
DCOLOR tempBins [preCompContrib -> nBins];
getPreCompContribByPhoton(&photon, photonIdx, preCompContrib,
tempBins
);
for (b = 0; b < preCompContrib -> nBins; b++) {
/* Accumulate binned contributions into rcontrib's LUT entry
for this modifier, weighted by ray coefficient. Also sum
total contributions. */
for (k = 0; k < 3; k++)
photonContrib [k] = tempBins [b] [k] *
rcData -> rayCoeff [k];
addcolor(modCont -> cbin [b], photonContrib);
addcolor(rcData -> totalContrib, photonContrib);
}
}
}
else {
/* NO BINNING; get single contribution directly from photon,
scale by ray coefficient and sum to total contribs */
getPhotonFlux(&photon, photonContrib);
multcolor(photonContrib, rcData -> rayCoeff);
addcolor(modCont -> cbin [0], photonContrib);
addcolor(rcData -> totalContrib, photonContrib);
}
}
return 0;
}
void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray,
COLOR totalContrib
)
/* Get precomputed light source contributions at ray -> rop for all
specified modifiers and accumulate them in pmapContribTab (which maps
to rcontrib's binning LUT). Also return total precomputed contribs. */
{
PreCompContribRCData rcData;
/* Ignore sources */
if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype))
return;
/* Get cumulative path coefficient up to photon lookup point */
raycontrib(rcData.rayCoeff, ray, PRIMARY);
setcolor(totalContrib, 0, 0, 0);
/* Rcontrib results passed to per-modifier child contrib pmaps */
rcData.ray = ray;
rcData.totalContrib = totalContrib;
/* Iterate over child contribution photon maps for each modifier and
* collect their contributions */
lu_doall(pmap -> preCompContribTab, getPreCompContribByMod, &rcData);
}
#endif /* PMAP_CONTRIB */
c4science · Help