Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmapdens.c
3772c0ebd22epmapmerge

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 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 1 / (4 * PI) */
scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
}
#ifdef PMAP_PHOTONFLOW
void lightFlowPhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Evaluate irradiance from volume photons at position 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 irrad) */
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);
/* Restore ray normal */
flipsurface(ray);
}
#endif /* PMAP_PHOTONFLOW */
c4science · Help