diff --git a/pmapdens.c b/pmapdens.c index 99a1aef..26b4666 100644 --- a/pmapdens.c +++ b/pmapdens.c @@ -1,251 +1,251 @@ #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" /* Photon map lookup functions per type */ void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = { photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity, photonDensity, photonDensity }; 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 +#ifdef PMAP_NONEFOUND sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", ray -> ro ? ray -> ro -> oname : "", 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); + getPhotonFlux(photon, flux); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on photon dist */ scalecolor(flux, 2 * (1 - sqn -> dist2 / r2)); -#endif +#endif addcolor(irrad, flux); } - + /* Divide by search area PI * r^2, 1 / PI required as ambient - normalisation factor */ + 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))); return; } #ifdef PMAP_PHOTONFLOW void volumePhotonSphIrrad (PhotonMap *pmap, RAY *ray, COLOR irrad) /* Evaluate 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". * The found photons in the search queue must originate from disparate * paths for an unbiased solution. This requires enabling the photon * path ID filter callback in findPhotons() to remove duplicate paths. */ { unsigned i; float r2; COLOR flux; Photon *photon; const PhotonSearchQueueNode *sqn; setcolor(irrad, 0, 0, 0); if (!pmap -> maxGather) return; findPhotons(pmap, ray); /* Issue warning if few too photons, as the photon path filter almost * invariably results in short lookups due to the reduced effective * photon density after dropping duplicates, thereby invalidating the * automatic search radius adjustment. As a workaround for now, it's * up to the user to the user to compensate this with a sufficiently * large fixed search radius. */ 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 */ scalecolor(flux, 2 * (1 - sqn -> dist2 / r2)); -#endif +#endif addcolor(irrad, flux); -#ifdef PMAP_DEBUGPATHS +#ifdef PMAP_DEBUGPATHS printf( - "%f\t%f\t%f\t%f\t%f\t%d:%010d\n", + "%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 +#endif } -#ifdef PMAP_DEBUGPATHS - printf("r = %f\n", sqrt(r2)); -#endif +#if defined(PMAP_DEBUGPATHS) + printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2)); +#endif /* Divide accumulated flux by sphere surface area and ambient normalisation factor PI */ scalecolor(irrad, 1 / (4 * PI * PI * r2)); return; } #endif