gsKdNode.h

Go to the documentation of this file.

# pragma once

#include <gsHSplines/gsAABB.h>

namespace gismo {

template<short_t d, class Z = index_t>
struct gsKdNode
{
    // Defines the type of the box
    typedef          gsAABB<d, Z> kdBox;
    typedef typename kdBox::point point;

    int axis; 

    Z pos;

    int level ;

    kdBox * box;

    gsKdNode * parent;

    gsKdNode * left;

    gsKdNode * right;

    gsKdNode() : axis(-2) ,level(0), box(0),
                 parent(0), left(0), right(0)
    { }

    gsKdNode(point const & upp) : axis(-1) , level(0),
                                  parent(0), left(0) , right(0)
    { 
        // Initial box, upp is expected to be indexed in finest level
        box = new kdBox( point::Zero(), upp);
    }

    gsKdNode(point const & low, point const & upp) : axis(-1) , level(0),
                                                     parent(0), left(0) , right(0)
    {
        // Initial box, upp is expected to be indexed in finest level
        box = new kdBox( low, upp);
    }

    gsKdNode(kdBox const & bb) : axis(-1) , level(0),
                                 parent(0), left(0) , right(0)
    { 
        // ..
    }
    
    gsKdNode(const gsKdNode & o, gsKdNode * parentNode = NULL) : axis(o.axis), level(o.level)
    {
        parent = parentNode;
        if ( axis == -1 )
        {
            GISMO_ASSERT( (o.left == 0) && (o.right == 0), 
                          "Problem: leaf with children." );
            box  = new kdBox(*o.box);
            left = right = NULL;
        }
        else
        {
            GISMO_ASSERT( o.box == 0, 
                          "Problem: split node with box." );
            pos   = o.pos;
            left  = new gsKdNode(*o.left , this);
            right = new gsKdNode(*o.right, this);
            box   = NULL;
        }
    }

    ~gsKdNode()
    {
        // TODO: non-recursive

        if ( isLeaf() ) 
        {
            delete box;
        }
        else
        {
            delete left;
            delete right;
        }
    }

    // Box Accessors
    const point & lowCorner() const 
    { 
        GISMO_ASSERT(box, "Asked for lowCorner at node without box data.");
        return box->first ; 
    }

    const point & uppCorner() const 
    { 
        GISMO_ASSERT(box, "Asked for uppCorner at node without box data.");
        return box->second; 
    }

    bool isLeaf() const { return axis == -1; }

    bool isRoot() const { return parent == NULL; }

    bool isTerminal() const 
    { return (axis!=-1) && (left->axis==-1) && (right->axis==-1); }

    bool isLeftChild()  const { return parent!=NULL && this==parent->left; }

    bool isRightChild() const { return parent!=NULL && this==parent->right; }

    bool isDegenerate() const
    { return (box->first.array() >= box->second.array()).any(); }

    gsKdNode * sibling() const
    { 
        GISMO_ASSERT( parent != 0, "Root does not have a sibling.");
        return (parent->left == this ? parent->right : parent->left ); 
    }

    void multiplyByTwo()
    {
        if ( isLeaf() )
        {
            box->first .array() *= 2;
            box->second.array() *= 2;
        }
        else
        {
            pos *= 2;
        }
    }

    void divideByTwo()
    {
        if ( isLeaf() )
        {
            box->first .array() /= 2;
            box->second.array() /= 2;
        }
        else
        {
            pos /= 2;
        }
    }

    inline void split()
    {
        GISMO_ASSERT( (left == 0) && (right == 0),
                      "Can only split leaf nodes.");
        GISMO_ASSERT( axis > -1, "Split axis not prescribed.");

        // Make new left and right children
        left          = new gsKdNode;
        right         = new gsKdNode;
        // Set axis to -1 (since they are leaves)
        left ->axis   =
        right->axis   = -1;
        // Set parent to this node
        left ->parent = 
        right->parent = this;
        // Set level
        left ->level  = 
        right->level  = level;
        // Set box
        left ->box    = box;    
        right->box    = new kdBox(*box);
        // Detach box from parent (is now at left child)
        box = NULL;
        // Resize properly the box coordinates
        left ->box->second[axis] = 
        right->box->first [axis] = pos;
    }

    inline void merge()
    {
        GISMO_ASSERT( (left->isLeaf()) && (right->isLeaf()),
                      "Can only merge terminal nodes.");

        // Recover box
        box = left->box;
        left->box = NULL;
        box->second[axis] = right->box->second[axis];
        axis  = - 1;
        level = left->level;

        // Delete children
        delete  left;
        left  = NULL;
        delete right;
        right = NULL;
    }


    void split(int splitAxis, Z splitPos)
    {
        GISMO_ASSERT( box->second[splitAxis] != splitPos, "Degenerate split " << box->second[splitAxis] <<" != "<<splitPos);
        GISMO_ASSERT( box->first [splitAxis] != splitPos, "Degenerate split " << box->first[splitAxis]  <<" != "<<splitPos);
        axis = splitAxis;
        pos  = splitPos;
        split();
    }

    // TODO: remove
    void nextMidSplit()
    {        
        axis = ( parent == 0 ? 0 : (parent->axis+1)%d );        
        pos  = box->first [axis] + 
            (box->second[axis] - box->first[axis])/2 ;
        split(); // Can be degenerate
    }

    void anyMidSplit(int index_level)
    {        
        const unsigned h = 1 << (index_level - level) ;
        const unsigned mask = ~(h - 1);
        for ( unsigned i = 0; i < d; ++i )
        {
            const unsigned c = 
                (box->first [i] + (box->second[i] - box->first[i])/2) & mask ;
            if ( c != box->first [i] ) // avoid degenerate split
            {
                split(i, c);
                return;
            }
        }
    }


    gsKdNode * adaptiveAlignedSplit(kdBox const & insBox, int index_level)
    {
        const unsigned h = 1 << (index_level - level) ;
        
        for (short_t i = 0; i < d; ++i)
        {
            const Z c1 = insBox. first[i] - insBox. first[i] % h; //floor
            const Z cc = insBox.second[i] % h;
            const Z c2 = insBox.second[i] + (cc ? h-cc : 0 ); // ceil

            if ( c1 > box->first[i] )
            {
                // right child intersects insBox
                split(i, c1 );
                return right;
            }
            else if ( c2 < box->second[i]  )
            {
                // left child intersects insBox
                split(i, c2 );
                return left;
            }
        }
        return NULL;
    }

    // TODO: remove
    gsKdNode * adaptiveSplit(kdBox const & insBox)
    {
        // assumption: insBox intersects box
        for ( unsigned i = 0; i < d; ++i )
        {
            // TODO: strategy: Try to split as close to the middle as possible.
            if ( insBox.first[i] > box->first[i] )
            {
                axis = i;
                pos  = insBox.first[i];
                split();
                return right;
            }
            else if ( insBox.second[i] < box->second[i] )
            {
                axis = i;
                pos  = insBox.second[i];
                split();
                return left ;
            }
        }
        // insBox is equal to this->box or they do not overlap
        return NULL;
    }

    friend std::ostream & operator<<(std::ostream & os, const gsKdNode & n)
    {
        if ( n.isLeaf() ) 
        {
            os << "Leaf node ("<< n.box->first.transpose() <<"), ("
               << n.box->second.transpose() <<"). level="<<n.level<<" \n";
        }
        else
        {
            os << "Split node, axis= "<< n.axis <<", pos="<< n.pos <<"\n";
        }

        return os;
    }

};


}// namespace gismo