src/SAHKDTree.cpp

Go to the documentation of this file.

#include <algorithm>
#include <map>

#include "SAHKDTree.h"
#include "Node.h"
#include "Event.h"
#include "SplitPlane.h"
#include "Primitive.h"


SAHKDTree::~SAHKDTree()
{
}


void SAHKDTree::buildTree(  const PrimitiveList &   list,
                            const Box &             bbox )
{
    KDTree::buildTree(list, bbox);

    // create sorted initial event list
    EventList events;
    PrimitiveList::const_iterator it = list.begin();
    for (; it != list.end(); ++it)
    {
        generateEvents(events, *it, bbox);
    }
    std::sort(events.begin(), events.end());

    // create the tree recursively
    buildTree(mRootNode, bbox, events, 0);

    LOG("Building kd-tree done");
}


void SAHKDTree::buildTree(  Node *          node,
                            const Box &     bbox,
                            EventList &     events,
                            unsigned int    depth ) const
{
    assert(node);
    assert(!node->left());
    assert(!node->right());

    // find best split plane
    float minCost = INFINITY;
    SplitPlane bestPlane = findBestPlane(minCost, node, bbox, events);

//    LOG("\n\ndepth: " << depth << " :: " << node->primitives_count() << " prims , " << events.size()
//            << " events \nbbox: [" << bbox.minVertex() << " , " << bbox.maxVertex() << "] , split dir/side/pos: "
//            << bestPlane.direction() << " " << bestPlane.side() << " " << bestPlane.position()
//            << "   cost: " << minCost);

    if (terminate(node, minCost, depth))
    {
        // make leaf node
        node->setSplitAxis(NO_AXIS);
        return;
    }
    assert(events.size() > 0);

    // setup node's split axis and plane
    node->setSplitPlane(bestPlane.position());
    node->setSplitAxis(bestPlane.direction());

    // create both children
    Node * left     = new Node();
    Node * right    = new Node();
    node->setLeftRight(left, right);

    // Event lists for children
    EventList eLeft;
    EventList eRight;

    // create event and primitive lists for children
    std::pair<Box, Box> boxes = classifyAndSplice(eLeft, eRight, node, bbox, events, bestPlane);

    // free some memory before recursing
    node->clear_primitives();
    events.clear();

    // setup subtrees
    buildTree(left, boxes.first, eLeft, depth + 1);
    buildTree(right, boxes.second, eRight, depth + 1);
}


bool SAHKDTree::terminate(  const Node *    current,
                            float           minCost,
                            unsigned int    depth ) const
{
    return KDTree::terminate(current, minCost, depth) || minCost > mKIntersect * current->primitives_count();
}


float SAHKDTree::SAH(   Side &              side,   // out param
                        const Box &         bbox,
                        const SplitPlane &  plane,
                        unsigned int        number_left,
                        unsigned int        number_right,
                        unsigned int        number_middle ) const
{
    float probability_left  = bbox.area(plane);
    float probability_right = 1.0f - probability_left;  // FIXME maybe calculate accurate probabilities
    float left_cost     = cost(probability_left, probability_right, number_left + number_middle, number_right);
    float right_cost    = cost(probability_left, probability_right, number_left, number_right + number_middle);

    if (left_cost < right_cost)
    {
        side = S_LEFT;
        return left_cost;
    }
    else
    {
        side = S_RIGHT;
        return right_cost;
    }
}



SplitPlane SAHKDTree::findBestPlane(    float &             minCost,    // out param
                                        const Node *        node,
                                        const Box &         bbox,
                                        const EventList &   events ) const
{
    // triangle counters
    unsigned int NL[3];
    unsigned int NP[3];
    unsigned int NR[3];

    // initialize counters
    for (int k=0; k<3; k++)
    {
        NL[k] = NP[k] = 0;
        NR[k] = node->primitives_count();
    }

    SplitPlane plane;
    minCost = INFINITY;

    // find best split plane out of all candidates
    EventList::const_iterator it = events.begin();
    while (it != events.end())
    {
        unsigned int p_end      = 0;
        unsigned int p_planar   = 0;
        unsigned int p_start    = 0;

        SplitPlane p = it->plane();

        // count different events untill the next plane is found
        while (it != events.end()
                && p.direction() == it->plane().direction()
                && p.position() <= it->plane().position()
                && it->type() == ET_END)
        {
            ++p_end;
            ++it;
        }
        while (it != events.end()
                && p.direction() == it->plane().direction()
                && p.position() <= it->plane().position()
                && it->type() == ET_PLANAR)
        {
            ++p_planar;
            ++it;
        }
        while (it != events.end()
                && p.direction() == it->plane().direction()
                && p.position() <= it->plane().position()
                && it->type() == ET_START)
        {
            ++p_start;
            ++it;
        }

        // incrementally update counters
        NP[p.direction()] = p_planar;
        NR[p.direction()] -= p_planar + p_end;

        // calculate costs
        Side side;
        float cost = SAH(side, bbox, p, NL[p.direction()], NR[p.direction()], NP[p.direction()]);

        // keep plane with minimal costs
        if (cost < minCost
                && bbox.minVertex()[p.direction()] + EPSILON < p.position()
                && p.position() < bbox.maxVertex()[p.direction()] - EPSILON)
        {
            plane       = p;
            minCost     = cost;
            plane.setSide(side);
        }

        // incrementally update counters
        NL[p.direction()] += p_planar + p_start;
        NP[p.direction()] = 0;
    }

    return plane;
}


void SAHKDTree::generateEvents( EventList &     events,     // out param
                                Primitive *     primitive,
                                const Box &     bbox ) const
{
    Box b = primitive->bounds().clip(bbox);
    for (int k=0; k<3; k++) // for all dimensions
    {
        // if bounding box is planar in current dimension
        if (b.minVertex()[k] >= b.maxVertex()[k] - EPSILON)
        {
            events.push_back(Event(primitive, SplitPlane(b.minVertex()[k], Axis(k)), ET_PLANAR));
        }
        else
        {
            events.push_back(Event(primitive, SplitPlane(b.minVertex()[k], Axis(k)), ET_START));
            events.push_back(Event(primitive, SplitPlane(b.maxVertex()[k], Axis(k)), ET_END));
        }
    }
}


std::pair<Box, Box> SAHKDTree::classifyAndSplice(   EventList &         eLeft,  // out param
                                                    EventList &         eRight, // out param
                                                    Node *              node,
                                                    const Box &         bbox,
                                                    const EventList &   events,
                                                    const SplitPlane &  plane ) const
{
    // create a map of primitives to markers for classification
    std::map<Primitive *, Side> prims;
    for (unsigned int i=0; i<node->primitives_count(); i++)
    {
        prims[node->primitive(i)] = S_BOTH;
    }
    assert(node->primitives_count() == prims.size());

    // generate event and primitive lists for children
    // old left-only and right-only events are sorted and stored before split iterator
    // newly generated both-side events are stored after split iterator

    // classify
    EventList::const_iterator it = events.begin();
    for (; it != events.end(); ++it)
    {
        if (it->type() == ET_END
                && it->plane().direction() == plane.direction()
                && it->plane().position() <= plane.position())
        {
            // left only
            prims[it->primitive()] = S_LEFT;
        }
        else if (it->type() == ET_START
                    && it->plane().direction() == plane.direction()
                    && it->plane().position() >= plane.position())
        {
            // right only
            prims[it->primitive()] = S_RIGHT;
        }
        else if (it->type() == ET_PLANAR
                    && it->plane().direction() == plane.direction())
        {
            if (it->plane().position() < plane.position()
                    || (it->plane().position() <= plane.position()
                        && plane.side() == S_LEFT))
            {
                // left only
                prims[it->primitive()] = S_LEFT;
            }
            if (it->plane().position() > plane.position()
                    || (it->plane().position() >= plane.position()
                        && plane.side() == S_RIGHT))
            {
                // right only
                prims[it->primitive()] = S_RIGHT;
            }
        }
        // both -> generate new events for clipped primitives (later)
    }

    // splice
    for (it = events.begin(); it != events.end(); ++it)
    {
        if (prims[it->primitive()] == S_LEFT)
        {
            eLeft.push_back(*it);
        }
        else if (prims[it->primitive()] == S_RIGHT)
        {
            eRight.push_back(*it);
        }
    }

    // save sizes of sorted parts
    unsigned int leftSize   = eLeft.size();
    unsigned int rightSize  = eRight.size();

    // bounding volumes for children FIXME do we need to calculate correct bounding volumes?
    std::pair<Box, Box> boxes = bbox.split(plane);

    // generate new events for overlapping primitives and copy primitives to children
    std::map<Primitive *, Side>::iterator mapit = prims.begin();
    for (; mapit != prims.end(); ++mapit)
    {
        if (mapit->second == S_LEFT)
        {
            node->left()->add(mapit->first);
        }
        else if (mapit->second == S_RIGHT)
        {
            node->right()->add(mapit->first);
        }
        else
        {
            generateEvents(eLeft, mapit->first, boxes.first);
            generateEvents(eRight, mapit->first, boxes.second);

            node->left()->add(mapit->first);
            node->right()->add(mapit->first);
        }
    }

//    LOG("event split: " << leftSize << " -- " << eLeft.size()-leftSize
//            << " --- " << eRight.size()-rightSize << " -- " << rightSize);

    // iterators between sorted and unsorted parts
    EventList::iterator leftBoth    = eLeft.begin() + leftSize;
    EventList::iterator rightBoth   = eRight.begin() + rightSize;

    // debug
    assert(eLeft.begin() <= leftBoth && leftBoth <= eLeft.end());
    assert(eRight.begin() <= rightBoth && rightBoth <= eRight.end());

    // sort new events and merge them with old events
    std::sort(leftBoth, eLeft.end());
    std::sort(rightBoth, eRight.end());
    std::inplace_merge(eLeft.begin(), leftBoth, eLeft.end());
    std::inplace_merge(eRight.begin(), rightBoth, eRight.end());

    return boxes;
}

---