Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmapdens.c
a8ee1830db88pmapmerge

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
, lightFlowDensity, transLightFlowDensity
#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 distance.
(see Eq. 4.1, p.76 in Silverman, "Density Estimation for
Statistics and Data Analysis", 1st Ed., 1986, and
Wann Jensen, "Realistic Image Synthesis using Photon Mapping") */
scalecolor(flux, 1 - sqn -> dist2 / r2);
#endif
addcolor(irrad, flux);
}
#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 */
scalecolor(irrad, 2 / (PI * PI * r2));
#else
/* Normalise accumulated flux by search area PI * r^2, including
RADIANCE-specific lambertian factor PI */
scalecolor(irrad, 1 / (PI * PI * r2));
#endif
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 4 * PI */
scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
}
#ifdef PMAP_PHOTONFLOW
void lightFlowDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Evaluate (static) irradiance from volume photons incident on plane
defined by normal ray -> ron at point ray -> rop, ignoring
participating medium. This evaluation interprets the volume photon
map as a representation of the "flow" of light, i.e. a physical light
field.
If sphericalIrrad == 0, the planar irradiance is evaluated with respect
to the plane defined by the ray normal ray -> ron.
If sphericalIrrad == 1, the mean spherical irradiance is evaluated,
and the normal ray -> ron is ignored. */
{
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 (has no effect if evaluating mean spherical irradiance) */
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 -> maxDist2Limit)
);
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 distance
(see Eq. 4.1, p.76 in Silverman, "Density Estimation for
Statistics and Data Analysis", 1st Ed., 1986) */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
if (!sphericalIrrad) {
/* Associate photons with plane defined by ray normal */
for (i = 0; i < 3; i++) {
/* Jitter photon direction since stored as 8-bit signed int */
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
}
#ifdef PMAP_EPANECHNIKOV
/* Normalise accumulated flux by Epanechnikov kernel integral in 3D
(see Eq. 4.4, p.76 in Silverman, "Density Estimation for Statistics
and Data Analysis", 1st Ed., 1986), include RADIANCE-specific
lambertian factor PI */
scalecolor(irrad, 15 / (8 * PI * PI * r2 * sqrt(r2)));
#else
/* Normalise accumulated flux by spherical volume to obtain irradiance,
include RADIANCE-specific irradiance factor PI */
scalecolor(irrad, 3 / (4 * PI * PI * r2 * sqrt(r2)));
#endif
if (sphericalIrrad)
/* Scale irradiance by 1/4; this corresponds to the "true" scalar
illuminance as defined in Eq. 13 in:
Mangkuto R.A., "Research note: The accuracy of the mean spherical
semi-cubic illuminance approach for determining scalar illuminance",
Lighting Research & Technology, 2020; 52: 151-158. */
scalecolor(irrad, 0.25);
#ifdef PMAP_PHOTONFLOW_HEMI
else
/* Hemispherical mode; lightflow photons which lie behind the
plane were filtered out (see filterPhoton()), so renormalise
for hemisphere (3/(2 PI)). This reduces bias near "hard"
boundaries (scene cube limits or actual geometry with opaque
surfaces), in which case no photons may be present behind
the plane. */
scalecolor(irrad, 2.0);
#endif
/* Restore ray normal */
flipsurface(ray);
}
void transLightFlowDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Evaluate transient irradiance from volume photons incident on plane
defined by normal ray -> ron at point ray->rop and time pmap->time,
ignoring participating medium. Additionally, the speed of light
(c) in units of the scene geometry is obtained from pmap->velocity.
Physically, this corresponds to finding the photons that would pass
through the lookup point at the given time.
If sphericalIrrad == 0, the planar irradiance is evaluated with
respect to the plane defined by the ray normal ray->ron.
If sphericalIrrad == 1, the mean spherical irradiance is evaluated,
and the normal ray->ron is ignored. */
{
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 (has no effect if evaluating mean spherical irradiance) */
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);
/* 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 distance
(see Eq. 4.1, p.76 in Silverman, "Density Estimation for
Statistics and Data Analysis", 1st Ed., 1986) */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
if (!sphericalIrrad) {
/* Associate photons with plane defined by ray normal */
for (i = 0; i < 3; i++) {
/* Jitter photon direction since stored as 8-bit signed int */
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
}
#ifdef PMAP_EPANECHNIKOV
/* Normalise accumulated flux by Epanechnikov kernel integral in 3D
(see Eq. 4.4, p.76 in Silverman, "Density Estimation for Statistics
and Data Analysis", 1st Ed., 1986), include RADIANCE-specific
lambertian factor PI */
scalecolor(irrad, 15 / (8 * PI * PI * r2 * sqrt(r2)));
#else
/* Normalise accumulated flux by spherical volume to obtain irradiance,
include RADIANCE-specific irradiance factor PI */
scalecolor(irrad, 3 / (4 * PI * PI * r2 * sqrt(r2)));
#endif
if (sphericalIrrad)
/* Scale irradiance by 1/4; this corresponds to the "true" scalar
illuminance as defined in Eq. 13 in:
Mangkuto R.A., "Research note: The accuracy of the mean spherical
semi-cubic illuminance approach for determining scalar illuminance",
Lighting Research & Technology, 2020; 52: 151-158. */
scalecolor(irrad, 0.25);
#ifdef PMAP_PHOTONFLOW_HEMI
else
/* Hemispherical mode; lightflow photons which lie behind the
plane were filtered out (see filterPhoton()), so renormalise
for hemisphere (3/(2 PI)). This reduces bias near "hard"
boundaries (scene cube limits or actual geometry with opaque
surfaces), in which case no photons may be present behind
the plane. */
scalecolor(irrad, 2.0);
#endif
/* Restore ray normal */
flipsurface(ray);
}
#endif /* PMAP_PHOTONFLOW */
c4science · Help