src/rcrt/SAHKDtree.hpp

Go to the documentation of this file.

#ifndef SAHKDTREE_HPP_
#define SAHKDTREE_HPP_

#include "Traceable.h"
#include <vector>
#include <stack>
#include <list>
#include <limits>
#include <iostream>



namespace rcrt
{

enum Side
{
        BOTH_SIDES = 0,
        LEFT_ONLY = 1,
        RIGHT_ONLY = 2,
};

enum OptimalCostSide
{
        NOT_SET = 0,
        LEFT = 1,
        RIGHT = 2,
};

enum EventType
{
        END = 0,
        IN_PLANE = 1,
        START = 2,
        
};

template <class T>
class KDNode {
public:
        KDNode<T>() {
                prim = 0;
                children[0] = 0;
                children[1] = 0;
        };
        
        ~KDNode() {
                if (children[0]) delete children[0];
                if (children[1]) delete children[1];
                if (prim) delete prim;
        }
        
        bool leaf, empty;
        KDNode<T>* children[2];
        Axis splitAxis;
        float plane;
        std::vector<T*>* prim;
        
};

template <class T> struct TraceableEntry;

template <class T>
struct PlanarEvent {
        TraceableEntry<T>* traceable;
        Axis axis;
        float plane;
        EventType type;
        OptimalCostSide optSide;
};

template <class T>
struct TraceableEntry {
        T* t;
        Side side;
};


template <class T>
class SAHKDtree : public rcrt::Traceable
{
private:
        KDNode<T>* root;
        AABB bbox;
        unsigned int maxDepth, leafSize;
public:
        SAHKDtree<T>() : root(0)
        {
                maxDepth = 40;
                leafSize = 2;
        }
        
        SAHKDtree<T>(std::vector<T*>* trList, int depth = 40, int lSize = 2) : root(0), maxDepth(depth), leafSize(lSize)
        {
                buildTree(trList);
        }
        
        virtual ~SAHKDtree() {
                
        }
        
        typedef typename std::vector< T* >::iterator vecTIter;
        typedef typename std::vector<TraceableEntry<T>* >::iterator vecTEIter;
        typedef typename std::list<PlanarEvent<T> >::iterator listPEIter;
        
        Intersection intersect(Ray& r) const
        {
                Intersection is = bbox.intersect(r);
                
                if (!is.isValid()) {
                        return is;
                }
                
                return intersectNode(root, r);          
        }
        
        Intersection intersectNode(KDNode<T>* n, Ray& r) const {
                
                Intersection is;
                
                if (n->empty) {
                        return is;
                }
                if (n->leaf) {
                        for (vecTIter i = n->prim->begin(); i != n->prim->end(); ++i) {
                                Intersection currentIs = (*i)->intersect(r);
                                if (currentIs.isValid() && currentIs.getDistance() < is.getDistance()) {
                                        is = currentIs;
                                }
                        }
                        return is;
                }
                
                float rayAxisDir = r.dir()[n->splitAxis];
                float rayAxisOrg = r.org()[n->splitAxis];
                float splitDist = (n->plane-rayAxisOrg)/rayAxisDir;
                
                int near = rayAxisDir < 0;
                
                if (r.minDist() - splitDist >  100 * std::numeric_limits<float>::epsilon() ) {
                        return intersectNode(n->children[!near],r);
                }
                
                if (splitDist - r.maxDist() > 100 * std::numeric_limits<float>::epsilon() ) {
                        return intersectNode(n->children[near],r);
                }
                
                float maxDist = r.maxDist();
                r.setMaxDist(std::min(splitDist + 100 * std::numeric_limits<float>::epsilon(),maxDist));
                Intersection nearIS = intersectNode(n->children[near],r);
                r.setMaxDist(maxDist);
                
                float minDist = r.minDist();
                r.setMinDist(std::max(splitDist - 100 * std::numeric_limits<float>::epsilon(),minDist));
                Intersection farIS = intersectNode(n->children[!near],r);
                r.setMinDist(minDist);
                
                if (nearIS.isValid() && nearIS.getDistance() < farIS.getDistance()) {
                        return nearIS;
                }
                else {
                        return farIS;
                }
                
        }
        
        const AABB& getBoundingBox() const
        {
                return bbox;
        }
        
        const Point3D& getCentroid() const
        {
                return bbox.getCentroid();
        }
        
        void setTraceables(std::vector<T*>* trList) {
                buildTree(trList);
        }
        
        void buildTree(std::vector<T*>* trList) {
                
                std::cout << "[SAHKDtree]: building tree size = "<< trList->size() << std::endl;
                
                if (root) delete root;
                
                bbox = calcBoundingBox(trList);
                
                
                std::vector<TraceableEntry<T>* > trEList;
                std::list<PlanarEvent<T> >* eventList = new std::list<PlanarEvent<T> >();
                
                for( vecTIter i = trList->begin(); i != trList->end(); ++i ) {
                        TraceableEntry<T>* trEntry = new TraceableEntry<T>();
                        trEntry->side = BOTH_SIDES;
                        trEntry->t = *i;
                        
                        trEList.push_back(trEntry);
                        createPlanarEvents(trEntry, eventList);
                }
                
                // sort event list
                eventList->sort(planarEventCompare);
                //std::cout << "Events Sorted" << std::endl;
                
                root = buildNode(bbox, trEList, eventList, 0);
                std::cout << "             built tree" << std::endl; 
        }
        
        KDNode<T>* buildNode(AABB& nbox, std::vector<TraceableEntry<T>* > trList, std::list<PlanarEvent<T> >* eventList, unsigned int depth) {
                
                //std::cout << "building Node - depth=" << depth << " size=" << trList.size() << " nbox=" << nbox << std::endl;
                
                // check wether to build a leaf
                if (trList.size() <= leafSize || depth >= maxDepth || nbox.getSurfaceArea() < 0.00001 || eventList->size() == 0) {
                        return buildLeaf(trList);
                }
                
                // determining best split
                
                int N = trList.size();
                
                int Nl[3], Np[3], Nr[3];
                
                for (int i = 0; i < 3; ++i) {
                        Nl[i] = 0; Np[i] = 0; Nr[i] = N;
                }
                
                //std::cout << "finding best split" << std::endl;
                
                PlanarEvent<T>* bestSplit = 0;
                float minCost = std::numeric_limits<float>::infinity();
                
                for (int k = 0; k < 3; ++k) {
                        listPEIter i = eventList->begin();
                        
                        while (i != eventList->end()) {
                                
                                // get next split candidate
                                
                                while (i != eventList->end() && (*i).axis != intToAxis(k)) {
                                        ++i;
                                }
                                
                                if (i == eventList->end()) break;
                                
                                PlanarEvent<T>* candidate = &(*i);
                                
                                // initializing p
                                int p[3] = {0};
                                
                                ++p[(*i).type];
                                
                                ++i;
                                
                                while (i != eventList->end() && (*i).plane == candidate->plane) {
                                        if ((*i).axis == intToAxis(k)) {
                                                ++p[(*i).type];
                                        }
                                        ++i;
                                }
                                
                                //moving plane on p
                                
                                //std::cout << "plane=" << p[IN_PLANE] << " start=" << p[START] << " end=" << p[END] << std::endl;
                                
                                Np[k] = p[IN_PLANE]; Nr[k] -= p[IN_PLANE]; Nr[k] -= p[END];
                                
                                float cost = SAH(candidate, nbox, Nl[k], Nr[k], Np[k]);
                                
                                //std::cout << "cost=" << cost << std::endl;
                                
                                if (cost < minCost) {
                                        minCost = cost;
                                        bestSplit = candidate;
                                }
                                
                                Nl[k] += p[START]; Nl[k] += p[IN_PLANE]; Np[k] = 0;
                        }
                }
                
                if (minCost >= 10 * trList.size()) {
                        return buildLeaf(trList);
                }
                
                //std::cout << "classification " << eventList->size() << std::endl;
                
                for (vecTEIter i = trList.begin(); i != trList.end(); ++i) {
                        (*i)->side = BOTH_SIDES;
                }
                
                //std::cout << "set both" << std::endl;
                
                for (listPEIter i = eventList->begin(); i != eventList->end(); ++i) {
                        //std::cout << "processing tri " << ((*i).axis) << " " << ((*i).traceable == 0) << " " << (bestSplit == 0) << std::endl;
                        if ((*i).type == END && (*i).axis == bestSplit->axis && (*i).plane <= bestSplit->plane) {
                                //std::cout << "LO1" << std::endl;
                                (*i).traceable->side = LEFT_ONLY;
                        }
                        else if ((*i).type == START && (*i).axis == bestSplit->axis && (*i).plane >= bestSplit->plane) {
                                //std::cout << "RO1" << std::endl;
                                (*i).traceable->side = RIGHT_ONLY;
                        }
                        else if ((*i).type == IN_PLANE && (*i).axis == bestSplit->axis) {
                                //std::cout << "IP" << std::endl;
                                if ((*i).plane < bestSplit->plane || ((*i).plane == bestSplit->plane && bestSplit->optSide == LEFT)) {
                                        //std::cout << "LO2" << std::endl;
                                        (*i).traceable->side = LEFT_ONLY;
                                }
                                if ((*i).plane > bestSplit->plane || ((*i).plane == bestSplit->plane && bestSplit->optSide == RIGHT)) {
                                        //std::cout << "RO2" << std::endl;
                                        (*i).traceable->side = RIGHT_ONLY;
                                }
                        }
                }
                
                //std::cout << "building tr lists" << std::endl;
                
                std::list<PlanarEvent<T> >* El = new std::list<PlanarEvent<T> >();
                std::list<PlanarEvent<T> >* Er = new std::list<PlanarEvent<T> >();
                std::vector<TraceableEntry<T>* > ltrs;
                std::vector<TraceableEntry<T>* > rtrs;
                
                for (vecTEIter i = trList.begin(); i != trList.end(); ++i) {
                        if ((*i)->side == LEFT_ONLY) {
                                ltrs.push_back(*i);
                        }
                        else if ((*i)->side == RIGHT_ONLY) {
                                rtrs.push_back(*i);
                        }
                        else {
                                ltrs.push_back(*i);
                                rtrs.push_back(*i);
                        }
                }
                
                //std::cout << "building exclusive lists" << std::endl;
                
                for (listPEIter i = eventList->begin(); i != eventList->end(); ++i) {
                        if ((*i).traceable->side == LEFT_ONLY) {
                                El->push_back(*i);
                        }
                        else if ((*i).traceable->side == RIGHT_ONLY) {
                                Er->push_back(*i);
                        }
                }
                
                // generating new Events for overlapping triangles
                
                //std::cout << "generate events for overlap" << std::endl;
                
                std::list<PlanarEvent<T> > Ebl;
                std::list<PlanarEvent<T> > Ebr;
                
                Point3D minP = nbox.getMinP();
                Point3D maxP = nbox.getMaxP();
                Point3D maxLP = maxP; maxLP[bestSplit->axis] = bestSplit->plane;
                Point3D minRP = minP; minRP[bestSplit->axis] = bestSplit->plane;
                
                AABB lBox(minP,maxLP);
                AABB rBox(minRP, maxP);
                
                //std::cout << "enter for 1" << std::endl;
                
                for (vecTEIter i = trList.begin(); i != trList.end(); ++i) {
                        if ((*i)->side == BOTH_SIDES) {
                                
//                              createPlanarEvents(*i, &Ebl, lBox);
//                              createPlanarEvents(*i, &Ebr, rBox);
//                              AABB lclipBox =(*i)->t->clipBox(lBox);
//                              AABB rclipBox =(*i)->t->clipBox(rBox);
//                              //std::cout << "lclipBox: " << lclipBox << " rclipBox: " << rclipBox << std::endl;
//                              if (!lclipBox.isEmpty()){
//                                      //std::cout << "lclipBox: " << lclipBox << std::endl;
//                                      createPlanarEvents(*i, &Ebl, lclipBox);
//                              }
//                              if (!rclipBox.isEmpty()) {
//                                      //std::cout << "rclipBox: " << rclipBox << std::endl;
//                                      createPlanarEvents(*i, &Ebr, rclipBox);
//                              }
                                AABB lclipBox;
                                AABB rclipBox;
                                //std::cout << "clip" << std::endl;
                                (*i)->t->clipPlane(bestSplit->axis, bestSplit->plane, lclipBox, rclipBox);
                                lclipBox.contract(lBox);
                                rclipBox.contract(rBox);
                                //std::cout << "lclipBox: " << lclipBox << " rclipBox: " << rclipBox << std::endl;
                                createPlanarEvents(*i, &Ebl, lclipBox);
                                //std::cout << "left events created" << std::endl;
                                createPlanarEvents(*i, &Ebr, rclipBox);
                                //std::cout << "right events created" << std::endl;
                        }
                }
                
                //std::cout << "sort and merge lists" << std::endl;
                
                // sorting lists of newly generated events
                
                Ebl.sort(planarEventCompare);
                Ebr.sort(planarEventCompare);
                
                // generate left and right event list
                
                El->merge(Ebl,planarEventCompare);
                Er->merge(Ebr,planarEventCompare);
                
                // create new node
                
                //std::cout << "Create Node" << std::endl;
                
                KDNode<T>* newNode = new KDNode<T>();
                newNode->plane = bestSplit->plane;
                newNode->splitAxis = bestSplit->axis;
                newNode->leaf = false;
                newNode->empty = false;
                newNode->children[0] = buildNode(lBox, ltrs, El, depth+1);
                newNode->children[1] = buildNode(rBox, rtrs, Er, depth+1);
                
                //std::cout << "done" << std::endl;
                
                delete eventList;
                
                return newNode;
        }
        
        KDNode<T>* buildLeaf(std::vector<TraceableEntry<T>* > trList) {
                
                //std::cout << "building leaf" << std::endl;
                
                KDNode<T>* leaf = new KDNode<T>();
                
                leaf->leaf =true;
                
                // check if leaf is empty
                if (trList.size() == 0) {
                        leaf->empty = true;
                        return leaf;
                }
                
                leaf->empty = false;
                
                leaf->prim = new std::vector<T*>();
                                
                // adding primitives to node
                for (vecTEIter i = trList.begin(); i != trList.end(); ++i) {
                        leaf->prim->push_back((*i)->t);
                }
                
                return leaf;
        }
        
        static bool planarEventCompare(PlanarEvent<T>& a, PlanarEvent<T>& b)
        {
                return a.plane < b.plane || (a.plane == b.plane && a.type < b.type);
        }
        
        
        void createPlanarEvents(TraceableEntry<T>* trEntry, std::list<PlanarEvent<T> >* eventList, AABB tbbox)
        {
                // creating planar events for every axis
                for (int axis = 0; axis < 3; ++axis) {
                        if (tbbox.getMin(intToAxis(axis)) == tbbox.getMax(intToAxis(axis))) {
                                // generating IN_PLANE event
                                PlanarEvent<T> planeEvent;
                                planeEvent.axis = intToAxis(axis);
                                planeEvent.plane = tbbox.getMin(intToAxis(axis));
                                planeEvent.traceable = trEntry;
                                planeEvent.type = IN_PLANE;
                                planeEvent.optSide = NOT_SET;
                                
                                eventList->push_back(planeEvent);
                        }
                        else {
                                // generating START and END events
                                PlanarEvent<T> startEvent;
                                startEvent.axis = intToAxis(axis);
                                startEvent.plane = tbbox.getMin(intToAxis(axis));
                                startEvent.traceable = trEntry;
                                startEvent.type = START;
                                startEvent.optSide = NOT_SET;
                                
                                eventList->push_back(startEvent);
                                
                                PlanarEvent<T> endEvent;
                                endEvent.axis = intToAxis(axis);
                                endEvent.plane = tbbox.getMax(intToAxis(axis));
                                endEvent.traceable = trEntry;
                                endEvent.type = END;
                                endEvent.optSide = NOT_SET;
                                
                                eventList->push_back(endEvent);
                        }
                }
        }
        
        void createPlanarEvents(TraceableEntry<T>* trEntry, std::list<PlanarEvent<T> >* eventList)
        {
                const AABB& tbbox = trEntry->t->getBoundingBox();
                createPlanarEvents(trEntry, eventList, tbbox);

        }
        
        float C(float& Pl, float& Pr, int Nl, int Nr)
        {       
                float tCost = 1;
                float iCost  = 10;
                float lambda = (Nl == 0 || Nr == 0) ? 0.8 : 1;
                return lambda * (tCost + iCost * (Pl * Nl + Pr * Nr));
        }
        
        float SAH(PlanarEvent<T>* p, AABB& V, int& Nl, int& Nr, int& Np) {
                Point3D minP = V.getMinP();
                Point3D maxP = V.getMaxP();
                Point3D maxLP = maxP; maxLP[p->axis] = p->plane;
                Point3D minRP = minP; minRP[p->axis] = p->plane;
                
                AABB Vl(minP,maxLP);
                AABB Vr(minRP, maxP);
                
                float vol = V.getSurfaceArea();
                
                //std::cout << "vol:" << vol << std::endl;
                
                float Pl = Vl.getSurfaceArea()/vol;
                float Pr = Vr.getSurfaceArea()/vol;
                
                float Cl = C(Pl,Pr,Nl+Np,Nr);
                float Cr = C(Pl,Pr,Nl,Np+Nr);
                
                //std::cout << "SAH - Cl=" << Cl << " Cr=" << Cr <<" Pl=" << Pl << " Pr=" << Pr << " Nl=" << Nl << " Nr=" << Nr << " Np=" << Np << std::endl;
                
                p->optSide = (Cl<Cr) ? LEFT : RIGHT;
                
                return (Cl<Cr) ? Cl : Cr;
        }
        
        static AABB calcBoundingBox(std::vector<T*>* trList)
        {
                AABB box;
                for(unsigned int i = 0; i < trList->size(); i++){
                        box.extend(trList->at(i)->getBoundingBox());
                }
                return box;
        }
        
        AABB clipBox(AABB& cbox) const {
                return getBoundingBox();
        }
        
        static Axis intToAxis(int i)
        {
                if (i==0) return X_AXIS;
                if (i==1) return Y_AXIS;
                if (i==2) return Z_AXIS;
                return X_AXIS;
        }
};

}

#endif /*SAHKDTREE_HPP_*/

---