Homec4science
Diffusion RADIANCE Photon Map (photonflow) Browse

RadiancePmap/pmapdens.c
39cad12f1bddphotonflow

Last Change
Raw File

pmapdens.c
View Options

#ifndef lint
static const char RCSid[] = "$Id$";
#endif
/*
======================================================================
Photon map density estimation routines
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Fraunhofer Institute for Solar Energy Systems,
supported by the German Research Foundation
(DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS))
(c) Lucerne University of Applied Sciences and Arts,
supported by the Swiss National Science Foundation
(SNSF #147053, "Daylight Redirecting Components")
(c) Tokyo University of Science,
supported by the JSPS Grants-in-Aid for Scientific Research
(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
======================================================================
$Id$
*/
#include "pmapdens.h"
#include "pmapbias.h"
#include "otypes.h"
#include "random.h"
/* Photon map lookup functions per type */
void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
photonDensity, photonDensity,
#ifdef PMAP_PHOTONFLOW
lightFlowPhotonDensity
#endif
};
void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Surface-bound photon density estimate for global and caustic photons.
Returns irradiance at ray -> rop. */
{
unsigned i;
float r2;
COLOR flux;
Photon *photon;
const PhotonSearchQueueNode *sqn;
setcolor(irrad, 0, 0, 0);
if (!pmap -> maxGather)
return;
/* Ignore sources */
if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype))
return;
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;
}
if (pmap -> minGather == pmap -> maxGather) {
/* No bias compensation. Just do a plain vanilla estimate */
sqn = pmap -> squeue.node + 1;
/* Average radius^2 between furthest two photons to improve accuracy */
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
/* Skip the extra photon */
for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
getPhotonFlux(photon, flux);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on photon dist */
scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
#endif
addcolor(irrad, flux);
}
/* Divide by search area PI * r^2, 1 / PI required as ambient
normalisation factor */
scalecolor(irrad, 1 / (PI * PI * r2));
return;
}
else
/* Apply bias compensation to density estimate */
biasComp(pmap, irrad);
}
void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
/* Surface-bound photon density estimate for (single) precomputed photon.
Returns precomputed irradiance at ray -> rop. */
{
Photon p;
setcolor(irrad, 0, 0, 0);
/* Ignore sources */
if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
return;
if (find1Photon(preCompPmap, r, &p))
/* p contains a found photon, so get its irradiance, otherwise it
* remains zero under the assumption all photons are too distant
* to contribute significantly */
getPhotonFlux(&p, irrad);
}
void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Volume photon density estimate using Henyey-Greenstein phase function.
Returns inscattered irradiance at ray -> rop. */
{
unsigned i;
float r2, gecc2, ph;
COLOR flux;
Photon *photon;
const PhotonSearchQueueNode *sqn;
setcolor(irrad, 0, 0, 0);
if (!pmap -> maxGather)
return;
findPhotons(pmap, ray);
/* Need at least 2 photons */
if (pmap -> squeue.tail < 2)
return;
gecc2 = ray -> gecc * ray -> gecc;
sqn = pmap -> squeue.node + 1;
/* Average radius^2 between furthest two photons to improve accuracy */
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
/* Skip the extra photon */
for (i = 1; i < pmap -> squeue.tail; i++, sqn++) {
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
/* Compute phase function for inscattering from photon */
if (gecc2 <= FTINY)
ph = 1;
else {
ph = DOT(ray -> rdir, photon -> norm) / 127;
ph = 1 + gecc2 - 2 * ray -> gecc * ph;
ph = (1 - gecc2) / (ph * sqrt(ph));
}
getPhotonFlux(photon, flux);
scalecolor(flux, ph);
addcolor(irrad, flux);
}
/* Divide by search volume 4 / 3 * PI * r^3 and phase function
normalization factor 1 / (4 * PI) */
scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
}
#ifdef PMAP_PHOTONFLOW
void volumePhotonSphIrrad (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Evaluate mean spherical irradiance from volume photons inscattered at
* ray->rop, ignoring participating medium. This evaluates the scalar
* irradiance of a physical light field represented by "photon flow". */
{
unsigned i;
float r2;
COLOR flux;
Photon *photon;
const PhotonSearchQueueNode *sqn;
setcolor(irrad, 0, 0, 0);
if (!pmap -> maxGather)
return;
findPhotons(pmap, ray);
/* Need at least 2 photons */
if (pmap -> squeue.tail < 2)
return;
/* Issue warning if few too photons; this will particularly happen
if PMAP_PHOTONFILT is enabled, since photons from duplicate paths
are pruned during the search, and this is not accounted for by
the automatic search radius adjustment. As a workaround for now,
it's up to the user provide a sufficiently large fixed search
radius with the -ma option to rtrace */
if (pmap -> squeue.tail < pmap -> squeue.len) {
sprintf(
errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
error(WARNING, errmsg);
}
sqn = pmap -> squeue.node + 1;
/* Average radius^2 between furthest two photons to improve accuracy */
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
/* Skip the extra photon */
for (i = 1; i < pmap -> squeue.tail; i++, sqn++) {
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
getPhotonFlux(photon, flux);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on photon dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(flux, 1 - sqn -> dist2 / r2);
#endif
addcolor(irrad, flux);
#ifdef PMAP_DEBUGPATHS
printf(
"%f\t%f\t%f\t%f\t%f\t%d:%010d\n",
sqrt(sqn -> dist2),
photon -> pos [0], photon -> pos [1], photon -> pos [2],
flux [0], photon -> proc, photon -> primary
);
#endif
}
#ifdef PMAP_DEBUGPATHS
printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2));
#endif
#ifdef PMAP_PATHFILT
/* Divide accumulated flux by sphere surface area and ambient
normalisation factor PI */
scalecolor(irrad, 1 / (4 * PI * PI * r2));
#else
/* Divide accumulated flux by sphere volume and ambient normalisation
factor PI */
scalecolor(irrad, 3 / (4 * PI * PI * r2 * sqrt(r2)));
#endif
}
void volumePhotonCubicIllum (PhotonMap *pmap, RAY *ray, COLOR illum)
/* Evaluate cubic illuminance according to Cuttle from volume photons
* inscattered at ray->rop, ignoring participating medium. The evaluation
* interprets the volume photon map as a representation of the physical
* light field ("photon flow"). */
{
unsigned ax, p, i ,j;
float r2, cosNorm, volNorm;
Photon *photon;
const PhotonSearchQueueNode *sqn;
COLOR photonFlux, E [6],
Evec [3], Esymm [3],
EvecLum, EsymmSum = {0, 0, 0};
FVECT norm0;
const FVECT norm [6] = {
/* Cube face normals (INVERTED since volume photon directions
point TOWARDS their halfspace of incidence */
{-1, 0, 0}, {1, 0, 0}, {0, -1, 0}, {0, 1, 0}, {0, 0, -1}, {0, 0, 1}
};
setcolor(illum, 0, 0, 0);
if (!pmap -> maxGather)
return;
/* Save ray normal */
VCOPY(norm0, ray -> ron);
#if 1
/* Perform separate lookup for each axis. This improves accuracy
particularly for low bandwitdhs, as the radius adjusts to the
density of incident photons along each axis. This comes at the
expense of 6 lookups.*/
for (ax = 0; ax < 6; ax++) {
/* Set cube face normal for current axis */
VCOPY(ray -> ron, norm [ax]);
setcolor(E [ax], 0, 0, 0);
/* Find photons incident from this axis, filtering by normal */
findPhotons(pmap, ray);
/* Need at least 2 photons, else skip */
if (pmap -> squeue.tail < 2)
continue;
if (pmap -> squeue.tail < pmap -> squeue.len) {
sprintf(
errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
error(WARNING, errmsg);
}
sqn = pmap -> squeue.node + 1;
/* Avg radius^2 between furthest two photons to improve accuracy */
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
#ifdef _PMAP_DEBUGPATHS
printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2));
#endif
/* Skip the extra photon */
for (p = 1; p < pmap -> squeue.tail; p++, sqn++) {
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
getPhotonFlux(photon, photonFlux);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
/* Project photon dir onto axis */
cosNorm = DOT(norm [ax], photon -> norm) / 127;
scalecolor(photonFlux, cosNorm);
addcolor(E [ax], photonFlux);
#ifdef PMAP_DEBUGPATHS
printf(
"%f\t%f\t%f\t%f\t%f\t%d:%010d\n", sqrt(sqn -> dist2),
photon -> pos [0], photon -> pos [1], photon -> pos [2],
photonFlux [0], photon -> proc, photon -> primary
);
#endif
}
/* Normalise illuminance for spherical volume */
scalecolor(E [ax], 3 / (4 * PI * r2 * sqrt(r2)));
printf(
"E [%s%d] = %.1f\n",
ax & 1 ? "-" : "+", ax >> 1, luminance(E [ax])
);
}
#else
/* Perform 1 lookup for all axes, and distribute photon flux
according incident direction. While faster than separate lookups
per axis, the accuracy is poor along those axes with few incident
photons, particularly in conjunction with low bandwidths. */
findPhotons(pmap, ray);
/* Need at least 2 photons */
if (pmap -> squeue.tail < 2)
return;
if (pmap -> squeue.tail < pmap -> squeue.len) {
sprintf(
errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
error(WARNING, errmsg);
}
sqn = pmap -> squeue.node + 1;
/* Avg radius^2 between furthest two photons to improve accuracy */
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
/* Volumetric normalisation factor */
#ifdef PMAP_PATHFILT
volNorm = 1 / (2 * PI * r2);
#else
# if 1
volNorm = 3 / (4 * PI * r2 * sqrt(r2));
#else
/* Hemispherical volume */
volNorm = 3 / (2 * PI * r2 * sqrt(r2));
#endif
#endif
#ifdef PMAP_DEBUGPATHS
printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2));
#endif
for (ax = 0; ax < 6; ax++) {
/* Set cube face normal for current axis */
setcolor(E [ax], 0, 0, 0);
}
/* Skip the extra photon */
for (p = 1; p < pmap -> squeue.tail; p++, sqn++) {
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
getPhotonFlux(photon, photonFlux);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
for (ax = 0; ax < 6; ax++) {
cosNorm = DOT(norm [ax], photon -> norm) / 127;
if (cosNorm > 0) {
/* Project photon dir onto axis */
scalecolor(photonFlux, cosNorm);
addcolor(E [ax], photonFlux);
}
}
#ifdef PMAP_DEBUGPATHS
printf(
"%f\t%f\t%f\t%f\t%f\t%d:%010d\n", sqrt(sqn -> dist2),
photon -> pos [0], photon -> pos [1], photon -> pos [2],
photonFlux [0], photon -> proc, photon -> primary
);
#endif
}
for (ax = 0; ax < 6; ax++) {
/* Normalise illuminance for spherical volume */
scalecolor(E [ax], volNorm);
printf(
"E [%s%d] = %.1f\n",
ax & 1 ? "-" : "+", ax >> 1, luminance(E [ax])
);
}
#endif
/* Calculate illuminance vector and symmetric illuminance. Unlike
Cuttle's original formulation, we defer scalar reduction and
preserve RGB components in chromatic analysis is required. */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
Evec [i][j] = E [i << 1][j] - E [(i << 1) + 1][j];
Esymm [i][j] = 0.5 * (
E [i << 1][j] + E [(i << 1) + 1][j] - fabs(Evec [i][j])
);
}
/* Reduce Evec to luminance per axis, sum Esymm over axes */
EvecLum [i] = luminance(Evec [i]);
addcolor(EsymmSum, Esymm [i]);
}
scalecolor(EsymmSum, 1.0 / 3);
printf(
"Evec\t= [%.1f\t%.1f\t%.1f]\n",
EvecLum [0], EvecLum [1], EvecLum [2]
);
printf(
"Esym\t= [%.1f\t%.1f\t%.1f]\n",
luminance(Esymm [0]), luminance(Esymm [1]), luminance(Esymm [2])
);
/* Reduce to scalars and return as triplet:
- illuminance vector magnitude |Evec|
- scalar symmetry |Esymm|
- scalar illuminance |Esr|. */
illum [0] = sqrt(DOT(EvecLum, EvecLum));
illum [1] = luminance(EsymmSum);
illum [2] = 0.25 * illum [0] + illum [1];
/* Normalise by RADIANCE-specific irradiance factor PI */
scalecolor(illum, 1 / PI);
/* Restore ray normal */
VCOPY(ray -> ron, norm0);
}
#endif /* PMAP_PHOTONFLOW */
#ifdef PMAP_PHOTONFLOW
void lightFlowPhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Evaluate irradiance from volume photons at ray->rop for plane defined
* by normal ray -> ron, ignoring participating medium. The evaluation
* interprets the volume photon map as a representation of the "flow"
* of light, i.e. a physical light field. */
{
unsigned p, i;
float r2, cosNorm;
Photon *photon;
const PhotonSearchQueueNode *sqn;
COLOR photonFlux;
FVECT photonDir;
setcolor(irrad, 0, 0, 0);
if (!pmap -> maxGather)
return;
/* Flip normal since volume photon directions point towards plane of
incidence */
flipsurface(ray);
/* Find photons with incident direction filtered by dot product with
normal (note that found photons may also be located _behind_ the
plane, as if having passed through it. */
findPhotons(pmap, ray);
if (pmap -> squeue.tail < pmap -> squeue.len) {
sprintf(errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
error(WARNING, errmsg);
}
/* Need at least 2 photons, else bail out */
if (pmap -> squeue.tail < 2)
return;
sqn = pmap -> squeue.node + 1;
/* Avg radius^2 between furthest two photons to improve accuracy */
r2 = max(sqn -> dist2, (sqn + 1) -> dist2);
r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2));
#ifdef PMAP_DEBUGPATHS
printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2));
#endif
/* Skip the extra photon */
for (p = 1; p < pmap -> squeue.tail; p++, sqn++) {
photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
getPhotonFlux(photon, photonFlux);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
/* Jitter photon direction since stored as 8-bit signed int */
for (i = 0; i < 3; i++) {
photonDir [i] = photon -> norm [i];
if (photon -> norm [i] && photon -> norm [i] & 127 < 127)
photonDir [i] += photon -> norm [i] & 0x80
? frandom() : -frandom();
}
/* Project photon direction onto normal */
cosNorm = DOT(ray -> ron, photonDir) / 127;
scalecolor(photonFlux, cosNorm);
addcolor(irrad, photonFlux);
#ifdef PMAP_DEBUGPATHS
printf("%f\t%f\t%f\t%f\t%f\t%d:%010d\n", sqrt(sqn -> dist2),
photon -> pos [0], photon -> pos [1], photon -> pos [2],
photonFlux [0], photon -> proc, photon -> primary
);
#endif
}
/* Normalise accumulated flux by spherical volume, include
RADIANCE-specific irradiance factor PI */
scalecolor(irrad, 3 / (4 * PI * PI * r2 * sqrt(r2)));
/* Restore ray normal */
flipsurface(ray);
}
#endif /* PMAP_PHOTONFLOW */
c4science · Help