PhotonMapping.cc

Go to the documentation of this file.

#include <esg/energy/PhotonMapEnergy.h>
#include <esg/explorer/NExtentsExplorer.h>
#include <gra/reflection/PhotonMapping.h>
#include <assert.h>

using namespace gra;

const int PhotonMapping::DIRECT_ILLUMINATION   = (1<<0);
const int PhotonMapping::SPECULAR_REFLECTION   = (1<<1);
const int PhotonMapping::CAUSTICS              = (1<<2);
const int PhotonMapping::INDIRECT_ILLUMINATION = (1<<3);
const int PhotonMapping::GLOBAL_ILLUMINATION   = (DIRECT_ILLUMINATION |
                                                  SPECULAR_REFLECTION |
                                                  CAUSTICS |
                                                  INDIRECT_ILLUMINATION);

const unsigned PhotonMapping::PARTITIONING      = 30;
const unsigned PhotonMapping::PARTITIONING_2    = PARTITIONING * PARTITIONING;
const double   PhotonMapping::ANTIBIAS_FRACTION = .001;

#if 0
bool PhotonMapping::_interpolate_irradiance(const vector<Photon*>& np,
                                            Color3f&               irradiance)
{
    double  sumW = 0.0;
    Color3f aux;
    const PhotonMap::Photon* p;
    Vector3 pos;
    double dist;
    double w;

    /*
     * E(x,n) = sum_[w_i>1/a](w_i(x,n) * E_i(x_i)) / sum_[w_i>1/a](w_i(x,n))
     */

    irradiance.set(0.0, 0.0, 0.0);
    
    for (int i = 0; i < np->size(); i++) {
        p = &((*np)[i]);
        if (!p->pIrradCache) continue;
        pos.set(np->pointLocation());
        pos.sub(p->location());
        dist = pos.length();
        if (p->pIrradCache->maxDistance <= dist) continue;
        w = p->pIrradCache->meanDistance / dist;
        sumW += w;
        aux.set(p->pIrradCache->irradiance);
        aux.scale(w);
        irradiance.add(aux);
    }

    if (sumW > 0.0) {
        irradiance.scale(1.0 / sumW);
        return true;
    } else
        return false;
}
#endif

bool PhotonMapping::_interpolate_irradiance(const vector<IrradianceCache::Value*>& np,
                                            Color3f&                irradiance)
{
    /*
     * E(x,n) = sum_[w_i>1/a](w_i(x,n) * E_i(x_i)) / sum_[w_i>1/a](w_i(x,n))
     */

    irradiance.set(0.0, 0.0, 0.0);
    double sumW = 0.0;
    
    for (unsigned i = 0; i < np.size(); i++) {
        float w = np[i]->getWeight();
        const Color3f& c = np[i]->getIrradiance();
        sumW += w;
        irradiance.x += c.x * w;
        irradiance.y += c.y * w;
        irradiance.z += c.z* w;
    }

    if (sumW > 0.0) {
        irradiance.scale(1.0 / sumW);
        return true;
    } else
        return false;
}

double* PhotonMapping::_power_square(const vector<Photon*>& np,
                                     const Vector3&         normal,
                                     BRDF&                  brdf,
                                     const MatVisitor&      matVisitor,
                                     double&                squareEnergy)
{
    squareEnergy = 0.0;

    double* square = new double [PARTITIONING_2];
    for (unsigned i = 0; i < PARTITIONING_2; i++) *(square+i) = 0.0;
    
    /*
     * Compute number of samples per each area
     */

    Vector2  uv;
    unsigned u,v;
    double   avgPower;
    double   cell = 1. / PARTITIONING;
    Vector3  pw;

    for (unsigned i = 0; i < np.size(); i++) {
        uv.set(brdf.dir2uv(matVisitor, normal, np[i]->getDirection()));

        // find coordinates of cell inside square
        u = (v = 0);
        for (unsigned j = 0; j < PARTITIONING; j++) {
            if (uv.x > cell) { uv.x -= cell; u++; }
            if (uv.y > cell) { uv.y -= cell; v++; }
        }

        pw.set(np[i]->getPower());
        avgPower = (pw.x + pw.y + pw.z) / 3.;
        *(square + (PARTITIONING * v) + u) += avgPower;
        squareEnergy += avgPower;
    }

    return square;
}

double PhotonMapping::_importance_sampling(const Vector3&    normal,
                                           BRDF&             brdf,
                                           const MatVisitor& matVisitor,
                                           double*           square,
                                           double            squareEnergy,
                                           Vector3&          nDir)
{
    /*
     * Build histogram from accumulated energies in cells
     * and select random cell proportionally to
     * the energies
     */
    
    div_t    qr;
    double   biasEllim = ANTIBIAS_FRACTION * squareEnergy;
    double   r1 = ESG_DBL_RAND;
    double   r  = r1 * squareEnergy; /* rescale from [0,1) to 
                                      * [0, total_energy_in_square) */
    double*  val = square;
    unsigned index = 0;
    do {
        if (*val == 0.) {
            /*
             * To avoid bias we alliminate all regions with
             * zero energy (probability) by adding a small
             * fraction of the overall energy stored within
             * the unique square to these regions.
             */
            squareEnergy += biasEllim;
            r += r1 * biasEllim; // scale
            *val = biasEllim;
        }
        if (r <= *val) break;
        r -= *(val++);
    } while (++index < PARTITIONING_2); // should never arrive
                
    qr = div((int)index, (int)PARTITIONING);

    /*
     * Get random numbers inside selected (random) cell
     */
    
    double s1 = (qr.quot + ESG_DBL_RAND) / PARTITIONING;
    double s2 = (qr.rem  + ESG_DBL_RAND) / PARTITIONING;

    /*
     * Get final random direction
     */
    
    brdf.importanceSample(matVisitor, normal, s1, s2, nDir, NULL);
    return (squareEnergy/(*(square+index)*PARTITIONING_2));
}

void PhotonMapping::_sample_diffuse(PointEnv&   surr,
                                    unsigned    actDepth,
                                    bool        firstObjHit,
                                    Color3f&    contribution,
                                    MatVisitor& matVisitor)
{
    //static int a = 0;
    //static int b = 0;
    
    /*
     * Precise multiple diffuse reflections using
     * Monte Carlo ray tracing with irradiance caching
     */

    // Sample diffuse BRDF
    BRDF * pBRDF = surr.pVisitableObj->diffuseBRDF();
    if (!pBRDF) pBRDF = _pDiffuseBRDF;
    if (!pBRDF) return;

    double avgDiff = matVisitor.avgDiffuse();

    if (avgDiff <= 0.0) return;

    /*
     * Try to interpolate duffuse contribution by inspecting irradiance cache
     */
    if (_pIrradCache) {
        vector<IrradianceCache::Value*> buffer;
        Color3f irradiance;
        try {
            _pIrradCache->getValue(surr.intersection,
                                   surr.normal,
                                   buffer);
            if (buffer.size() > 0 && _interpolate_irradiance(buffer,irradiance)) {
                //cerr << "Vyuziti cache (interpolace/vypocet): " << ((float)(++a)/(float)b) << endl;
                contribution.add(irradiance);
                return;
            }
        } catch (out_of_range& ore) {
            cerr << ore.what() << endl;
        }
    }

    //cerr << "Vyuziti cache (interpolace/vypocet): " << ((float)a/(float)(++b)) << endl;
    
    /*
     * Irradiance caching failed.
     * Compute precise multiple diffuse reflections using Monte Carlo
     * ray tracing.
     *
     * Cast new ray(s) to get new irradiance value and possibly store
     * irradiance in cache
     */
    
    Vector3   nDir;
    Color3f   contrib(0.0, 0.0, 0.0);
    unsigned  castPhotons = ((_pIrradCache)
                             ? PARTITIONING_2
                             : _photonsPerSurface);
    Vector3   diff(matVisitor.diffuse());
    Vector3   refl;
    PointEnv* pNewSurr;
    double    meanDist = 0.0;
    unsigned  meanDistRays = 0;
    double*   lightingSquare = NULL;
    double    squareEnergy;

    if (_incomingLighting) {
        /*
         * Get closest photons from global photon map. These photons are
         * used for computation of average incoming lighting.
         */
        vector<Photon*>* np = _pGlobalMap->getNeighbourhood(surr.intersection,
                                                            &(surr.normal),
                                                            _maxDistToEst,
                                                            _numPhotonsToEst);
        if (np) {
            /*
             * Computer incoming energy distribution partitioned
             * in quare area
             */
            lightingSquare = _power_square(*np, surr.normal, *pBRDF,
                                           matVisitor, squareEnergy);
            delete np;
        }
    }

    if (!lightingSquare) {
        /*
         * Either incoming lighting is off or we got no photons from
         * global map. In both cases the scale factor is constant for
         * all casted rays.
         */
        refl.set(diff);            // scale irradiance by BRDF reflectance
        refl.scale(1.0 / avgDiff); // scale due to spectral differences
    }

    //cerr << _photonsPerSurface << endl;
    for (unsigned i = 0; i < castPhotons; i++) {
        if (ESG_DBL_RAND >= avgDiff) continue; // Russian roulette - absorption

        if (lightingSquare) {
            /*
             * Importance sampling.
             * Sample rays with respect to incoming lighting.
             *
             * => get the scale factor which is specific for each concrete
             *    photon and then get new random direction with respect to
             *    incoming lighting
             */
            refl.set(diff);
            refl.scale(_importance_sampling(surr.normal, *pBRDF, matVisitor,
                                            lightingSquare, squareEnergy,
                                            nDir));
        } else {
            /*
             * Get new random direction with respect to only shape of BRDF
             *
             * To do: stratified sampling (?)
             */
            pBRDF->importanceSample(matVisitor, surr.normal,
                                    ESG_DBL_RAND, ESG_DBL_RAND,
                                    nDir, NULL);
        }
        
        /*
         * Compute irradiance by casting the ray
         */
        if ((pNewSurr = _cast_ray(surr.intersection, nDir))) {
            assert(pNewSurr->mask & ENV_HAVE_DISTANCE);
            pNewSurr->energy.set(surr.energy);
            pNewSurr->energy.x *= refl.x;
            pNewSurr->energy.y *= refl.y;
            pNewSurr->energy.z *= refl.z;
            pNewSurr->mask |= ENV_HAVE_ENERGY;
            if (((pNewSurr->energy.x +
                  pNewSurr->energy.y +
                  pNewSurr->energy.z) / 3.) > _minRayWeight) {
                _illuminate(*pNewSurr, actDepth+1, firstObjHit, contrib);
            }

            if (_pIrradCache) {
                /*
                 * Compute harmonic mean distance in addition to irradiance
                 */     
                if (pNewSurr->distance != 0.0) {
                    meanDist += 1.0 / pNewSurr->distance;
                    meanDistRays ++;
                }
            }
            delete pNewSurr;
        } else {
            /*
             * Ray run away the scene. To append the backgound color as
             * contribution or to leave the contribution untached?
             * There is the question.
             *
             * Adding color of background seem to be correct at first sight.
             * But might produce color artifacts on edges of polygonal
             * surfaces. This is becaus the contribution is not attenuated
             * anyhow.
             */
            //contrib.add(_background);
        }

    } // for

    contrib.scale(1. / castPhotons);

    if (lightingSquare) delete [] lightingSquare;

    if (_pIrradCache) {
        /*
         * Store precisely computed irradiance to cache
         */
        if (meanDist > 0.0) {
            /*
             * Store irradiance in the closest photon
             */
            meanDist = meanDistRays / meanDist; 
            //meanDist = (1.0 / meanDistRays) / meanDist;
            try {
                _pIrradCache->addValue(surr.intersection,
                                       surr.normal,
                                       contrib,
                                       meanDist);
                //cout << _pIrradCache->getNumStoredValues() << endl;
            } catch(out_of_range& e) {
                cerr << e.what() << endl;
            }
            //contribution.set(1.0, 1.0, 1.0); // TEST !!!
        } else {
            /*
             * Wow, infinite harmonic mean distance. In this case we cannot
             * cache the irradiance but we can use it as the iradiance of
             * the point of interest.
             */
            
            //contribution.set(0.0, 1.0, 0.0); // TEST !!!
        }

    } //else contribution.set(0.0, 0.0, 0.0); // TEST !!!

    contribution.add(contrib);
}

void PhotonMapping::_sample_specular(PointEnv&   surr,
                                     unsigned    actDepth,
                                     bool        firstObjHit,
                                     Color3f&    contribution,
                                     MatVisitor& matVisitor)
{
    Vector3   nDir;
    Color3f   contrib(0.0, 0.0, 0.0);
    Vector3   refl;
    BRDF    * pBRDF = NULL;

    /*
     * We computes a (precise) specular and glossy reflections using
     * the Monte Carlo ray tracing
     */

    // => Sample specular BRDF
    pBRDF = surr.pVisitableObj->specularBRDF();
    if (!pBRDF) pBRDF = _pSpecularBRDF;
    if (!pBRDF) return;

    /*
     * Do russian roulette for every photon 
     */
    for (unsigned i = 0; i < _photonsPerSpecularSurface; i++) {
        double r1 = ESG_DBL_RAND;
        double rr = r1; // make a copy

        if ((rr -= matVisitor.avgSpecular()) < 0.) { // sample BRDF
            Vector3 prDir;
            _reflection_dir(surr, prDir); // prefered dir.
            pBRDF->importanceSample(matVisitor, prDir, r1, ESG_DBL_RAND,
                                    nDir, NULL);
            refl.set(matVisitor.specular());
            refl.scale(1. / matVisitor.avgSpecular());
            _cast_the_ray(surr, actDepth, true, contrib, nDir, refl);
        } else if ((rr -= matVisitor.avgReflection()) < 0.) { // reflection
            _reflection_dir(surr, nDir);
            refl.set(matVisitor.reflection());
            refl.scale(1. / matVisitor.avgReflection());
            _cast_the_ray(surr, actDepth, true, contrib, nDir, refl);
        } else if ((rr -= matVisitor.avgTransparency()) < 0.) { // refraction
            bool foh = ((_refraction_dir(surr,
                                         matVisitor,
                                         firstObjHit,
                                         nDir)
                         == RayTracing::REFR_DIR_REFLECTION)
                        ? false
                        : !firstObjHit);
            refl.set(matVisitor.transparency());
            refl.scale(1. / matVisitor.avgTransparency());
            _cast_the_ray(surr, actDepth, foh, contrib, nDir, refl);
        } // else absorption
    }
    
    contrib.scale(1./_photonsPerSpecularSurface);
    contribution.add(contrib);
}

void PhotonMapping::_illuminate(PointEnv& surr,
                                unsigned  actDepth,
                                bool      firstObjHit,
                                Color3f&  contribution)
{
    /*
     * We suppose that the given "surr" contains proper point of intersection,
     * normal vector and associated primitive
     *
     * Transform point of intersection and properly orient normal
     */
    _check_point_of_intersection(surr);

    /*
     * Inspect attributes that are necessary for local illumination,
     * reflection and transmission. Set correct BRDF and emittance for local
     * reflection. If the point of interest is the second point of refraction
     * ray then previous material is used (the object was not changed)
     */
    MatVisitor* pMatVisitor = _check_material(surr.pVisitableObj, firstObjHit);

    /*
     * Inspect photon maps
     */
    AutoPtr<EnergyCoat>* pAMap = surr.pVisitableObj->getEnergyState();
    if (pAMap) {
        EnergyCoat* pCoat = pAMap->origObject();
        if (pCoat) {
            PhotonMapEnergy* pE = dynamic_cast<PhotonMapEnergy*>(pCoat);
            if (pE) {
                _pGlobalMap = pE->getGlobalMap();
                _pCausticsMap = pE->getCausticsMap();
                if (_cacheIrradiance) {
                    _pIrradCache = pE->getIrradianceCache();
                    //cerr << _sceneDistance << endl;
                    if (!_pIrradCache && _sceneDistance > 0) {
                        //cerr << _maxError << endl;
                        if (pE->buildIrradianceCache(1.0f/_maxError,_sceneDistance))
                            _pIrradCache = pE->getIrradianceCache();
                        else
                            cerr << "PhotonMapping::_illuminate(): Irradiance cache will not be used." << endl;
                    }
                } else
                    _pIrradCache = NULL;
            }
        }
    }

    /*
     * If depth is zero then we always want the precise computation
     * We also want to check the illumination strategy on the only this
     * first level of ray traversal. In the next ray iterations
     * the illumination strategy is not checked and a full global illumination
     * is used instead.
     */
    if (actDepth == 0) { // always precise
        // 1) Direct illumination
        if (_illumStage & DIRECT_ILLUMINATION)
            _direct_illumination(surr, *pMatVisitor, actDepth, contribution);
        
        // 2) Caustics -- from caustics photon map
        if (_pCausticsMap && (_illumStage & CAUSTICS)) {
            Vector3 irrad;
            if (_pCausticsMap->getIrradiance(surr.intersection,
                                             NULL, /*&(surr.normal)*/
                                             _maxDistToEst,
                                             _numPhotonsToEst,
                                             irrad))
            {
                contribution.add(irrad);
            }
        }
        if (_illumStage & INDIRECT_ILLUMINATION) {
            // 3) Multiple diffuse reflections -- Monte Carlo ray tracing
            bool d = _diffReflection;
            _diffReflection = true;
            _sample_diffuse(surr, actDepth, firstObjHit,
                            contribution, *pMatVisitor);
            _diffReflection = d;
        }

        if (_illumStage & SPECULAR_REFLECTION) {
            // 4) Specular and glossy reflecion -- Monte Carlo ray tracing
            _specReflections++;
            _sample_specular(surr, actDepth, firstObjHit,
                             contribution, *pMatVisitor);
            _specReflections--;
        }

        //if (_pIrradCache) _pIrradCache->__debug(cerr);

        return;
    }

    /*
     * Check precision of computation
     */
    bool precise = false;
    if (_diffReflection) { // at least one diffuse reflection in path
        if ((surr.mask & ENV_HAVE_DISTANCE) &&
            (surr.distance < _preciseIllumDist))
            precise = true; // but close enought
    } else {  // no diffuse surface in path
        if (_specReflections <= _preciseIllumRefl) {
            precise = true; // either only a few specular reflections
        } else if ((surr.mask & ENV_HAVE_DISTANCE) &&
                   (surr.distance < _preciseIllumDist))
            precise = true; // or close enought
    }
    // to do: approximate if the ray contributes little to the pixel radiance

    /*
     * Direct illumination
     */
    if (precise) {
        // 1) Direct illumination
        _direct_illumination(surr, *pMatVisitor, actDepth, contribution);
            
        // 2) Caustics -- from caustics photon map
        if (_pCausticsMap) {
            Vector3 irrad;
            if (_pCausticsMap->getIrradiance(surr.intersection,
                                             &(surr.normal),
                                             _maxDistToEst,
                                             _numPhotonsToEst,
                                             irrad))
            {
                contribution.add(irrad);
            }
        }
        if (actDepth < _maxDepth) {
            // 3) Multiple diffuse reflections -- Monte Carlo ray tracing
            bool d = _diffReflection;
            _diffReflection = true;
            _sample_diffuse(surr, actDepth, firstObjHit,
                            contribution, *pMatVisitor);
            _diffReflection = d;

            // 4) Specular and glossy reflecion -- Monte Carlo ray tracing
            _specReflections++;
            _sample_specular(surr, actDepth, firstObjHit,
                             contribution, *pMatVisitor);
            _specReflections--;
        }
    } else {
        // 1) Direct illumination -- from global photon map if exists
        if (_pGlobalMap) {
            Vector3 irrad;
            if (_pGlobalMap->getIrradiance(surr.intersection,
                                           &(surr.normal),
                                           _maxDistToEst,
                                           _numPhotonsToEst,
                                           irrad))
            {
                contribution.add(irrad);
            }
        } else
            _direct_illumination(surr,*pMatVisitor,actDepth,contribution);
        
        // 2) Caustics -- from the global photon map above

        if (_diffReflection) {
            // 3) Multiple diffuse reflections
            //       -- energy from global map + stop recursion
        } else {
            if (actDepth < _maxDepth) {
                // 4) Specular and glossy reflection -- Monte Carlo ray tracing
                _specReflections++;
                _sample_specular(surr, actDepth, firstObjHit,
                                 contribution, *pMatVisitor);
                _specReflections--;
            }
        }
    }
}



// --------- public ----------------

Color3f* PhotonMapping::illuminatePoint(PointEnv & env)
{
    _specReflections = 0;
    _diffReflection  = false;
    _pGlobalMap      = NULL;
    _pCausticsMap    = NULL;
    _pIrradCache     = NULL;
    return PathTracing::illuminatePoint(env);
}

void PhotonMapping::setScene(Scene * s)
{
    PathTracing::setScene(s);

    float a[3][3] = {{1.,0.,0.}, {0.,1.,0.}, {0.,0.,1.}};

    if (_cacheIrradiance && s && s->root()) {
        float minD = FLT_MAX;
        float maxD = FLT_MIN;
        NExtentsExplorer e(a, 3);
        e.explore(*(s->root()));
        for (int i = 0; i < 3; i++) {
            Interval ext = e.result(i);
            if (ext.max > maxD) maxD = ext.max;
            if (ext.min < minD) minD = ext.min;
        }
        _sceneDistance = maxD - minD;
    }
}

---