src/rcrt/gi/PKDNode.cpp

Go to the documentation of this file.

#include "PKDNode.h"

using namespace std;

namespace rcrt{

PKDNode::PKDNode(std::list<PKDSplitPlane*>* splits, const float& maxRadius, const AABB& bbox )
{
        if(splits->size() == 3){
                //generate leaf, there is only one photon left with 3 planes
                //cout << "generating leaf" << endl;
                children[0] = 0;
                children[1] = 0;
                photon = splits->back()->p;
                splitAxis = splits->back()->axis;
                //finally delete the no longer needed splitplanes
                for(list<PKDSplitPlane* >::iterator iter = splits->begin(); iter != splits->end(); ++iter){
                        delete (*iter);
                }               
                delete splits;
                return;
        }
        //getting the best split photon requires evaluation of the cost
        //function C(p) = P(V_l)*N_l + P(V_r)*N_r, where P(V_i) = Vol(V_i)/Vol(V)
        //As Vol(V) is irrelevant for comparison of C(p1) < C(p2) we can
        //safely ignore it, thus P(V_i) = Vol(V_i) => C(p) = V_l * N_l
        //We seperate into C_l(p) = V_l*N_l and C_r(p) = V_r*N_r for
        //calculation.
        const int N = splits->size() / 3;//the list is always a multiple of 3
        int N_l[3] = {0,0,0};//number of photons to the left, separate for each axis
        int N_r[3] = {N-1,N-1,N-1};//          ''            right(minus the current)
        const Point3D& bMin(bbox.getMinP());
        const Point3D& bMax(bbox.getMaxP());
        PKDSplitPlane* bestSplit = 0;//track the best split
        float Cmin = numeric_limits<float>::infinity();
        
        for(list<PKDSplitPlane* >::const_iterator iter = splits->begin(); iter != splits->end(); ++iter){
                PKDSplitPlane* split = *iter;
                const Axis splitAxis = split->axis;
                const float& pos = split->p->getPos()[splitAxis];
                
                float cost = 0;
                
                if(N_l[splitAxis] != 0){
                        //approximating the Volume
                        float V_l = 1;
                        for(int k = 0; k < 3; k++){
                                if(k == splitAxis)
                                        //the volume extends along splitAxis to pos in the max
                                        V_l *= (pos-bMin[splitAxis] + 2 * maxRadius);
                                else
                                        V_l *= (bMax[k]-bMin[k] + 2 * maxRadius);
                        }
                        cost += V_l * N_l[splitAxis];
                }
                
                if(N_r[splitAxis] != 0){
                        //approximating the Volume
                        float V_r = 1;
                        for(int k = 0; k < 3; k++){
                                if(k == splitAxis)
                                        //the volume extends along splitAxis to pos in the min
                                        V_r *= (bMax[splitAxis] - pos + 2 * maxRadius);
                                else
                                        V_r *= (bMax[k]-bMin[k] + 2 * maxRadius);
                        }
                        cost += V_r * N_r[splitAxis];                   
                }
                //update best split and cost if this one is better
                if(cost < Cmin){
                        Cmin = cost;
                        bestSplit = split;
                }
                
                //update N_l and N_r
                N_l[splitAxis]++;
                N_r[splitAxis]--;
        }
        
        
        std::list<PKDSplitPlane* >* leftList = new std::list<PKDSplitPlane* >();
        std::list<PKDSplitPlane* >* rightList = new std::list<PKDSplitPlane* >();
        
        // now that we got us our split plane initialize this node
        photon = bestSplit->p;
        splitAxis = bestSplit->axis;
        
        // generate new bounding boxes that don't include this photon
        AABB leftBox;
        AABB rightBox;
        
        // and iterate once more over the list to seperate it into two
        const float splitValue = photon->getPos()[splitAxis];   
        for(list<PKDSplitPlane* >::const_iterator iter = splits->begin(); iter != splits->end(); ++iter){
                PKDSplitPlane* split = *iter;
                //do not include split planes of the photon in this node
                if(split->p == photon)
                        continue;
                
                //compare with split on split axis and put in correct list
                if(split->p->getPos()[splitAxis] < splitValue){
                        leftList->push_back(split);
                        //TODO can do the aabb more efficiently..3/2*N see PKDTree
                        leftBox.extend(split->p->getPos());
                }else{
                        rightList->push_back(split);
                        rightBox.extend(split->p->getPos());
                }
        }
        //delete unneeded list of splits for this node
        delete splits;
        
        // construct children
        if(leftList->size() != 0)
                children[0] = new PKDNode(leftList, maxRadius, leftBox);
        else
                children[0] = 0;
        
        if(rightList->size() != 0)
                children[1] = new PKDNode(rightList, maxRadius, rightBox);
        else
                children[1] = 0;
}

PKDNode::~PKDNode()
{
        if(children[0]) delete children[0];
        if(children[1]) delete children[1];
}

bool PKDNode::isLeaf() const
{
        return children[0] == 0 && children[1] == 0;
}

void PKDNode::getKNearest(std::vector<PKDSearchP>* result, const unsigned int& k,
                                const Point3D& pos, const Vec3D& normal, const float& maxRadiusSqr,
                                const float& flat) const
{
        const Point3D& pPos(photon->getPos());
        const float dx(pPos[0]-pos[0]);
        const float dy(pPos[1]-pos[1]);
        const float dz(pPos[2]-pos[2]);
        float distSqr(dx*dx+dy*dy+dz*dz);
        if(flat != 0){
                const float dotP(normal * (pPos - pos));
                const float ellipse(dotP * distSqr * 16 * flat);
                distSqr += ellipse;
        }
        
        float newMaxRadiusSqr = maxRadiusSqr;
        if(distSqr < maxRadiusSqr && normal * photon->getDir() < 0){
                // we found a photon that we want to insert
                //TODO only take it into account if it is inside an ellipsoid?
                if(result->size() == k-1){
                        //need to make it into a heap now
                        result->push_back(PKDSearchP(photon, distSqr));
                        make_heap(result->begin(), result->end());
                        //size is now k, retaining heap property from now on
                        //updating maxRadius, if k Photons found
                        newMaxRadiusSqr = result->front().distSqr;
                } 
                else if(result->size() >= k){
                        //need to retain heap property here
                        //check distance to photon at the front of the queue
                        if(distSqr < result->front().distSqr) {
                                //we need to put this photon in, pop the other
                                pop_heap(result->begin(), result->end());
                                result->pop_back();
                                result->push_back(PKDSearchP(photon, distSqr));
                                push_heap(result->begin(), result->end());
                        }
                        newMaxRadiusSqr = result->front().distSqr;
                } else {
                        // well, lets get it in there
                        result->push_back(PKDSearchP(photon, distSqr));
                        //push_heap(result->begin(),result->end());//for lazy, uncomment
                }
                
        }
        
        if(!isLeaf()){
                const float dist1D = pPos[splitAxis]-pos[splitAxis];
                int near = (dist1D > 0) ? 0 : 1; //which child is nearer?
                //search near subtree first
                if(children[near] != 0)
                        children[near]->getKNearest(result, k, pos, normal, newMaxRadiusSqr, flat);
                //search far subtree only if inside our search radius
                if(dist1D*dist1D <= newMaxRadiusSqr)
                        if(children[1-near] != 0)
                                children[1-near]->getKNearest(result, k, pos, normal, newMaxRadiusSqr, flat);
        }
}

}

---