src/rcrt/gi/PhotonMap.cpp

Go to the documentation of this file.

#include "PhotonMap.h"
#include <boost/thread/thread.hpp>
#include <boost/bind.hpp>
#include "../Ray.h"
#include "../materials/Material.h"

using namespace std;

namespace rcrt
{

PhotonMap::PhotonMap(Scene* s, const float& maxRad, const int& maxB, const float& flatten):
        scene(s),globalTree(maxRad * maxRad), causticTree(maxRad * maxRad),
        maxRadius(maxRad), maxBounces(maxB), currentRadius(maxRad / 8.0f),flat(flatten)
{
}

PhotonMap::~PhotonMap()
{
        for(unsigned int i = 0; i < photons.size(); i++)
                delete photons[i];
}

void PhotonMap::generate(const int& noPhotons)
{
        photons.reserve(noPhotons);
        caustics.reserve(noPhotons/4);
        vector<Light*>* lights = scene->getLights();
        //compute how many photons to send from each light
        //send more from brighter lights
        int noPhotonsPL[lights->size()];
        float totalPower = 0;
        int phots = 0;
        for(unsigned int i = 0; i < lights->size(); i++){
                totalPower += lights->at(i)->getPower().avg();
        }
        for(unsigned int i = 0; i < lights->size(); i++){
                noPhotonsPL[i] = (int)floor((lights->at(i)->getPower().avg() / totalPower) * noPhotons);
                phots += noPhotonsPL[i];
        }
        phots = noPhotons - phots;
        unsigned int i = 0;
        while(phots > 0){//distribute remaining photons
                if(i >= lights->size())
                        i = 0;
                noPhotonsPL[i]++;
                phots--;
        }
        
        RGBColor totalPowerRem(0);
        
        //multi-threaded sending of photons 
        i = 0;
        //compute first light extra if there is an uneven number of lights
        if(lights->size() % 2 != 0){
                vector<Photon*> caustics1;
                caustics1.reserve(noPhotonsPL[0]/4);
                vector<Photon*> photons1;
                photons1.reserve(noPhotonsPL[0]);
                int storedTotal = 0;
                cout << "[PhotonMap] Sending "<<noPhotonsPL[0]<<" photons for light: " << 0 << endl;
                boost::thread first(boost::bind(boost::mem_fn(&PhotonMap::generateJob), this, lights->at(0),
                                &photons1, &caustics1, boost::cref(noPhotonsPL[0]), boost::ref(storedTotal)));  
                i++;
                first.join();
                const float weight = 1.0f / float(storedTotal);
                cout << "[PhotonMap]    Sent " << noPhotonsPL[0] << " photons, with weight " << weight << endl;
                for(unsigned int k = 0; k < caustics1.size(); k++){
                        Photon* p = caustics1[k];
                        p->setWeight(weight);
                        caustics.push_back(p);
                        totalPowerRem = totalPowerRem + p->getPower();
                }
                for(unsigned int k = 0; k < photons1.size(); k++){
                        Photon* p = photons1[k];
                        p->setWeight(weight);
                        photons.push_back(photons1[k]);
                        //cout << "[PhotonMap] photon " << k << " p=" << p->getPower() << " b=" << p->getBounces() << endl;
                        totalPowerRem = totalPowerRem + p->getPower();
                }
        }
        //work two threads & light sources at the same time
        for(; i < lights->size(); i += 2){
                cout << "[PhotonMap] Sending "<< noPhotonsPL[i] << " photons for light: " << i << endl;
                cout << "[PhotonMap] Sending "<< noPhotonsPL[i+1] << " photons for light: " << i+1 << endl;
                vector<Photon*> caustics1;//reserve some space for caustic photons
                caustics1.reserve(noPhotonsPL[i]/4);
                vector<Photon*> caustics2;
                caustics2.reserve(noPhotonsPL[i+1]/4);
                vector<Photon*> photons1;
                photons1.reserve(noPhotonsPL[i]);
                vector<Photon*> photons2;
                photons2.reserve(noPhotonsPL[i+1]);
                int storedTotal1 = 0;
                int storedTotal2 = 0;
                
                boost::thread first(boost::bind(boost::mem_fn(&PhotonMap::generateJob), this, lights->at(i),
                                &photons1, &caustics1, boost::cref(noPhotonsPL[i]), boost::ref(storedTotal1)));
                boost::thread second(boost::bind(boost::mem_fn(&PhotonMap::generateJob), this, lights->at(i+1),
                                &photons2, &caustics2, boost::cref(noPhotonsPL[i+1]), boost::ref(storedTotal2)));
                
                first.join();
                const float weight1 = 1.0f / float(storedTotal1);
                cout << "[PhotonMap]    Sent " << noPhotonsPL[i] << " photons with weight "<< weight1 << endl;
                second.join();
                const float weight2 = 1.0f / float(storedTotal2);
                cout << "[PhotonMap]    Sent " << noPhotonsPL[i+1] << " photons with weight "<< weight2 << endl;
                for(unsigned int k = 0; k < caustics1.size(); k++){
                        Photon* p = caustics1[k];
                        p->setWeight(weight1);
                        caustics.push_back(p);
                        totalPowerRem = totalPowerRem + p->getPower();
                }
                for(unsigned int k = 0; k < caustics2.size(); k++){
                        Photon* p = caustics2[k];
                        p->setWeight(weight2);
                        caustics.push_back(p);
                        totalPowerRem = totalPowerRem + p->getPower();
                }
                for(unsigned int k = 0; k < photons1.size(); k++){
                        Photon* p = photons1[k];
                        p->setWeight(weight1);
                        photons.push_back(p);
                        totalPowerRem = totalPowerRem + p->getPower();
                }
                for(unsigned int k = 0; k < photons2.size(); k++){
                        Photon* p = photons2[k];
                        p->setWeight(weight2);
                        photons.push_back(p);
                        totalPowerRem = totalPowerRem + p->getPower();
                }
        }
        //got all the photons that we want, generate the trees
        cout << "          building global photon kdtree of size "<< photons.size() << endl;
        globalTree.buildTree(&photons);
        cout << "          building caustic photon kdtree of size " << caustics.size() << endl;
        causticTree.buildTree(&caustics);
        cout << "     total power in the photon map: " << totalPowerRem << endl;
}

void PhotonMap::generateJob(Light* light, vector<Photon*>* photons, vector<Photon*>* caustics,
                const int& noPhotons, int& storedTotal)
{
        int storedGlobal = 0;
        int storedCaustic = 0;
        while(storedGlobal+storedCaustic < noPhotons){
                //cout << " sending number "<< i << endl;
                Photon* photon = new Photon();
                //get photon from light
                light->emitPhoton(photon);
                //cout << "  emitted photon " << photon->getPos() << photon->getDir() << endl;
                //trace photon
                Ray photRay(photon->getPos(), photon->getDir());
                Intersection is(scene->intersect(photRay));
                //cout << "  traced photon ray" << endl;
                if(!is.isValid()){
                        //does not hit anything
                        //get a new photon
                        delete photon;
                        continue;
                } else {
                        complex<float> prevIOR(1,0);
                        is.setLastIOR(prevIOR);
                        Material* prevMat = is.getPrimitive()->getMaterial();
                        RGBColor prevPower = photon->getPower();
                        Vec3D prevIncDir = photon->getDir();
                        bool prevBackSide = is.backSide();
                        ScatterEvent scEvent = prevMat->scatterPhoton(is, photon);
                        ScatterEvent lastScEvent = scEvent;
                        if(scEvent == INVALID || scEvent == STORE){
                                delete photon;
                                continue;
                        }
                        //photon->addBounce();
                        bool causticPath = false;
                        
                        while(scEvent != INVALID
                                        && scEvent != STORE){
                                if(photon->getBounces() == maxBounces){
                                        scEvent = STORE;
                                        break;
                                }
                                //store a photon at the current intersection
                                if(photon->getBounces() > 2 && scEvent == DIFFUSE) {
                                        if(causticPath){
                                                caustics->push_back(
                                                        new Photon(prevPower,
                                                                        prevIncDir, photon->getPos(), photon->getBounces()));
                                                storedCaustic++;                                                
                                        } else {
                                                photons->push_back(
                                                        new Photon(prevPower,
                                                                        prevIncDir, photon->getPos(), photon->getBounces()));
                                                
                                                storedGlobal++;
                                        }
                                }
                                
                                causticPath = causticPath || scEvent == SPECULAR || scEvent == REFRACTED;
                                
                                Ray nextRay = Ray(photon->getPos(), photon->getDir());
                                is = scene->intersect(photRay);
                                if(!is.isValid()){//terminate path
                                        scEvent = INVALID;
                                        break;
                                }
                                /*if(prevMat->refracts() && !prevBackSide){
                                        //we traveled through a refracting medium to this point
                                        //attenuate the power
                                        const RGBColor& absorb(prevMat->getAbsorbance());
                                        const float& dist(is.getDistance());
                                        photon->scalePower(RGBColor(exp(-dist * absorb.r()),
                                                        exp(-dist * absorb.g()),
                                                        exp(-dist * absorb.b())));
                                }
                                prevBackSide = is.backSide();*/
                                if(lastScEvent == REFRACTED)
                                        is.setLastIOR(prevMat->getIOR());
                                else
                                        is.setLastIOR(complex<float>(1,0));
                                prevMat = is.getPrimitive()->getMaterial();
                                photRay = nextRay;
                                prevPower = photon->getPower();
                                prevIncDir = photon->getDir();
                                lastScEvent = scEvent;
                                scEvent = prevMat->scatterPhoton(is, photon);
                                photon->addBounce();
                                //don't follow photons that are too dark
                                if(photon->getPower().sum()-0.001f <= numeric_limits<float>::epsilon()){
                                        scEvent = INVALID;
                                        break;
                                }
                        }
                        if(scEvent == INVALID){
                                // we hit something useless / the path has terminated
                                // without getting absorbed or 
                                // store with previous values (incoming)
                                photon->setDir(prevIncDir);
                                photon->setPower(prevPower);
                                //store the photon
                                photons->push_back(photon);
                                storedGlobal++;
                        } else {
                                photons->push_back(photon);
                                storedGlobal++;
                        }                       
                }
        }
        
        storedTotal = storedGlobal + storedCaustic;
}

RGBColor PhotonMap::getRadiance(const Vec3D& wOut,
                const int& noGPhotons, const int& noCPhotons, Intersection& is)
{       
        Material* mat = is.getPrimitive()->getMaterial();
        Vec3D normal(mat->getShadingNormal(is));
        const Point3D& pos(is.getPosition());   
        
        //float currRSqr = currentRadius*currentRadius;
        float currRSqr = maxRadius*maxRadius;
        //bool tooFew = false;
        std::vector<PKDSearchP> global;
        if(noGPhotons != 0) globalTree.getKNearest(&global, noGPhotons, pos, normal, currRSqr, flat);
        
        std::vector<PKDSearchP> caustic;
        if(noCPhotons != 0) causticTree.getKNearest(&caustic, noCPhotons, pos, normal, currRSqr/2, flat);
        
        static const float gauss_norm = 1.8790;
        static const float alpha = 0.918 * gauss_norm;
        static const float beta = 1.953;
        static const float divG = 1.0f / (1.0 - exp(-beta));
        
        //cout << "hurz" << endl;
        float maxGR = -numeric_limits<float>::infinity();
        
        for(unsigned int i = 0; i < global.size(); i++){
                const float& distSqr(global[i].distSqr);
                //cout << "distSqr " << distSqr << endl;
                if(distSqr > maxGR) maxGR = distSqr;            
        }
        
        RGBColor globalR(0);
        for(unsigned int i = 0; i < global.size(); i++){
                const float& distSqr(global[i].distSqr);
                const float gauss = alpha * (1.0 - (1.0 - exp( -beta * distSqr / (M_PI * maxGR) )) * divG);
                Photon* p = global[i].p;
                globalR = globalR + p->getPower() * fabs(normal * p->getDir())
                                        * p->getWeight() * mat->sample(wOut, p->getDir(), is)
                                        * gauss;
        }
        
        float maxCR = -numeric_limits<float>::infinity();
        
        for(unsigned int i = 0; i < caustic.size(); i++){
                const float& distSqr(caustic[i].distSqr);
                if(distSqr > maxCR) maxCR = distSqr;
        }
        
        RGBColor causticR(0);
        for(unsigned int i = 0; i < caustic.size(); i++){
                const float& distSqr(caustic[i].distSqr);
                const float gauss = alpha * (1.0 - (1.0 - exp( -beta * distSqr / (M_PI * maxCR) )) * divG);
                Photon* p = caustic[i].p;
                causticR = causticR + p->getPower() * fabs(normal * p->getDir())
                                        * p->getWeight() * mat->sample(wOut, p->getDir(), is)
                                        * gauss;
        }
        //cout << globalR << causticR << endl;
        return globalR  + causticR;
}

const float& PhotonMap::getMaxRadius() const
{
        return maxRadius;
}

}

---