Page MenuHomec4science
Files F63768784

pmapcontrib.c
No OneTemporary
Actions

File Metadata

Created
Wed, May 22, 08:58

pmapcontrib.c
View Options

#ifndef lint
static const char RCSid[] = "$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $";
#endif
/*
=========================================================================
Photon map routines for precomputed light source contributions; these
routines interface to mkpmap.
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Lucerne University of Applied Sciences and Arts,
supported by the Swiss National Science Foundation
(SNSF #147053, "Daylight Redirecting Components",
SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment")
=========================================================================
$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $
*/
#include "pmapcontrib.h"
#include "pmapdiag.h"
#include "pmaprand.h"
#include "pmapmat.h"
#include "pmapsrc.h"
#include "pmapio.h"
#include "pmutil.h"
#include "otspecial.h"
#include "otypes.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 & (1 << 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 = 1 << 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, char *binVal, 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 (!binVal)
/* Fall back to single bin if no value specified */
binVal = "0";
eBinVal = eparse(binVal);
if (eBinVal -> type == NUM) {
/* Bin value is constant */
binCnt = (int)(evalue(eBinVal) + 1.5);
if (binCnt != 1) {
sprintf(errmsg, "invalid non-zero constant for bin %s", binVal);
error(USER, errmsg);
}
}
else if (binCnt <= 0) {
sprintf(errmsg, "undefined/invalid bin count for modifier %s", mod);
error(USER, errmsg);
}
/* Allocate and init contributions */
contrib = (MODCONT *)malloc(
sizeof(MODCONT) + sizeof(DCOLOR) * (binCnt - 1)
);
if (!contrib)
error(SYSTEM, "out of memory in addContribModifier()");
contrib -> outspec = NULL;
contrib -> modname = mod;
contrib -> params = binParm;
contrib -> nbins = binCnt;
contrib -> binv = eBinVal;
contrib -> bin0 = 0;
memset(contrib -> cbin, 0, sizeof(DCOLOR) * binCnt);
lutEntry -> data = lutEntry -> data = (char *)contrib;
++(*numContribs);
return(contrib);
}
void addContribModfile (LUTAB *contribTab, unsigned *numContribs,
char *modFile, char *binParm, char *binVal, int binCnt
)
/* Add light source modifiers from file modFile to contribution lookup table
* contribTab, and update numContribs.
* NOTE: This code is adapted from rcontrib.c */
{
char *mods [MAXMODLIST];
int i;
/* Find file and store strings */
i = wordfile(mods, MAXMODLIST, getpath(modFile, getrlibpath(), R_OK));
if (i < 0) {
sprintf(errmsg, "can't open modifier file %s", modFile);
error(SYSTEM, errmsg);
}
if (*numContribs + i >= MAXMODLIST - 1) {
sprintf(errmsg, "too many modifiers (max. %d) in file %s",
MAXMODLIST - 1, modFile
);
error(INTERNAL, errmsg);
}
for (i = 0; mods [i]; i++)
/* Add each modifier */
addContribModifier(contribTab, numContribs, mods [i],
binParm, binVal, binCnt
);
}
static int contribSourceBin (LUTAB *contribs, const RAY *ray)
/* Map contribution source ray to its bin for light source ray -> rsrc,
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)
/* We're 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 contribution source
* spawned photons (i.e. has srcIdx >= 0), it's appended to contribSrcHeap.
* If contribSrcRay == NULL, the current contribution source is still
* flushed, but no new source is set. Returns updated contribution source
* counter pmap -> numContribSrc. */
{
if (!pmap || !contribSrcHeap)
return 0;
/* Check if last contribution source has spawned photons (srcIdx >= 0,
* see newPhoton()), in which case we save it to the heap file before
* clobbering it. (Note this is short-term storage, so we doan' need
* da portable I/O stuff here). */
if (pmap -> lastContribSrc.srcIdx >= 0) {
if (!fwrite(&pmap -> lastContribSrc,
sizeof(PhotonContribSource), 1, contribSrcHeap
))
error(SYSTEM,
"failed writing photon contrib source in newPhotonContribSource()"
);
pmap -> numContribSrc++;
if (!pmap -> numContribSrc ||
pmap -> numContribSrc > PMAP_MAXCONTRIBSRC
);
}
/* Mark this contribution source unused with a negative source index
until its path spawns a photon (see newPhoton()) */
pmap -> lastContribSrc.srcIdx = -1;
/* Map ray to bin in anticipation that this contribution source will
be used, since the ray will be lost once a photon is spawned */
pmap -> lastContribSrc.srcBin = contribSourceBin(
pmap -> contribTab, contribSrcRay
);
return pmap -> numContribSrc;
}
PhotonContribSourceIdx buildContribSources (PhotonMap *pmap,
FILE **contribSrcHeap, char **contribSrcHeapFname,
PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps
)
/* Consolidate per-subprocess contribution sources heaps into array
* pmap -> contribSrc. Returns offset for contribution source index
* linearisation in pmap -> numContribSrc. The heap files in contribSrcHeap
* are closed on return. */
{
PhotonContribSourceIdx heapLen;
unsigned heap;
if (!pmap || !contribSrcHeap || !contribSrcOfs || !numHeaps)
return 0;
pmap -> numContribSrc = 0;
for (heap = 0; heap < numHeaps; heap++) {
contribSrcOfs [heap] = pmap -> numContribSrc;
if (fseek(contribSrcHeap [heap], 0, SEEK_END) < 0)
error(SYSTEM,
"failed photon contribution source seek in buildContribSources()"
);
pmap -> numContribSrc += heapLen = ftell(contribSrcHeap [heap]) /
sizeof(PhotonContribSource);
if (!(pmap -> contribSrc = realloc(pmap -> contribSrc,
pmap -> numContribSrc * sizeof(PhotonContribSource)
)))
error(SYSTEM,
"failed photon contrib source alloc in buildContribSources()"
);
rewind(contribSrcHeap [heap]);
if (fread(pmap -> contribSrc + contribSrcOfs [heap],
sizeof(PhotonContribSource), heapLen, contribSrcHeap [heap]
) != heapLen)
error(SYSTEM,
"failed reading photon contrib source in buildContribSources()"
);
fclose(contribSrcHeap [heap]);
unlink(contribSrcHeapFname [heap]);
}
return pmap -> numContribSrc;
}
static MODCONT *photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Locate photons near ray -> rop which originated from the light source
modifier indexed by ray -> rsrc, and accumulate their contributions
in pmap -> srcContrib. Return total contributions in irrad and pointer
to binned contributions (or NULL if none were found). */
{
const OBJREC *srcMod = findmaterial(source [ray -> rsrc].so);
MODCONT *contrib = (MODCONT*)lu_find(
pmap -> contribTab, srcMod -> oname
) -> data;
Photon *photon;
COLOR flux;
PhotonSearchQueueNode *sqn;
float r2, invArea;
unsigned i;
/* Zero bins for this source modifier */
memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins);
setcolor(irrad, 0, 0, 0);
if (!contrib || !pmap -> maxGather)
/* Modifier not in LUT or zero bandwidth */
return NULL;
/* Lookup photons */
pmap -> squeue.tail = 0;
/* Pass light source index to filter in findPhotons() via
pmap -> lastContribSrc */
pmap -> lastContribSrc.srcIdx = ray -> rsrc;
findPhotons(pmap, ray);
/* Need at least 2 photons */
if (pmap -> squeue.tail < 2) {
#ifdef PMAP_NONEFOUND
sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
ray -> ro ? ray -> ro -> oname : "<null>",
ray -> rop [0], ray -> rop [1], ray -> rop [2]
);
error(WARNING, errmsg);
#endif
return NULL;
}
/* Average radius^2 between furthest two photons to improve accuracy
and get inverse search area 1 / (PI * r^2). */
/* XXX: OMIT extra normalisation factor 1/PI for ambient calculation? */
sqn = pmap -> squeue.node + 1;
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
invArea = 1 / (PI * PI * r2);
/* Skip the extra photon */
for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
/* Get photon's contribution to density estimate */
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
getPhotonFlux(photon, flux);
scalecolor(flux, invArea);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on photon distance */
scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
#endif
addcolor(irrad, flux);
addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux);
}
return contrib;
}
void preComputeContrib (PhotonMap *pmap)
/* Precompute contributions for a random subset of finalGather *
pmap -> numPhotons photons, and build the precomputed contribution photon
map, discarding the original photons. */
/* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */
{
unsigned long i, numPreComp;
unsigned j;
PhotonIdx pIdx;
Photon photon;
RAY ray;
MODCONT *binnedContribs;
PhotonMap nuPmap;
repComplete = numPreComp = finalGather * pmap -> numPhotons;
if (verbose) {
sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n",
numPreComp
);
eputs(errmsg);
#if NIX
fflush(stderr);
#endif
}
/* Copy photon map for precomputed photons */
memcpy(&nuPmap, pmap, sizeof(PhotonMap));
/* Zero counters, init new heap and extents */
nuPmap.numPhotons = 0;
initPhotonHeap(&nuPmap);
for (j = 0; j < 3; j++) {
nuPmap.minPos [j] = FHUGE;
nuPmap.maxPos [j] = -FHUGE;
}
/* Record start time, baby */
repStartTime = time(NULL);
#ifdef SIGCONT
signal(SIGCONT, pmapPreCompReport);
#endif
repProgress = 0;
photonRay(NULL, &ray, PRIMARY, NULL);
ray.ro = NULL;
for (i = 0; i < numPreComp; i++) {
do {
/* Get random photon from stratified distribution in source heap
* to avoid duplicates and clustering */
pIdx = firstPhoton(pmap) + (unsigned long)(
(i + pmapRandom(pmap -> randState)) / finalGather
);
getPhoton(pmap, pIdx, &photon);
/* Init dummy photon ray with intersection and normal at photon
position. Set emitting light source index from origin. */
VCOPY(ray.rop, photon.pos);
ray.rsrc = photonSrcIdx(pmap, &photon);
for (j = 0; j < 3; j++)
ray.ron [j] = photon.norm [j] / 127.0;
/* Get contributions at photon position; retry with another photon
* if no contribs found */
binnedContribs = photonContrib(pmap, &ray, ray.rcol);
} while (!binnedContribs);
/* Append photon to new heap from ray */
newPhoton(&nuPmap, &ray);
/* TODO: Apply wavelet xform to binned contribs, apply lossy compression,
encode to RGBE */
#if 0
encodeContrib(binnedContribs, pmap -> compRatio);
#endif
/* Update progress */
repProgress++;
if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
pmapPreCompReport();
#ifdef SIGCONT
else signal(SIGCONT, pmapPreCompReport);
#endif
}
/* Flush heap */
flushPhotonHeap(&nuPmap);
#ifdef SIGCONT
signal(SIGCONT, SIG_DFL);
#endif
/* Trash original pmap, replace with precomputed one */
deletePhotons(pmap);
memcpy(pmap, &nuPmap, sizeof(PhotonMap));
if (verbose) {
eputs("\nRebuilding precomputed contribution photon map\n");
#if NIX
fflush(stderr);
#endif
}
/* Rebuild underlying data structure, destroying heap */
buildPhotonMap(pmap, NULL, NULL, 1);
/* TODO: Save wavelet coeffs here? */
}
#else
#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>\n", argv [0]);
fputs("Missing 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 = 1 << (l << 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);
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

Event Timeline

c4science · Help