BSPTreePoly.cpp

Go to the documentation of this file.

#include "moab/CartVect.hpp"
#include "moab/BSPTreePoly.hpp"
#include <cassert>
#include <cstdlib>
#include <set>

#undef DEBUG_IDS

namespace moab
{

struct BSPTreePoly::Vertex : public CartVect
{
    Vertex( const CartVect& v )
        : CartVect( v ), usePtr( 0 ), markVal( 0 )
#ifdef DEBUG_IDS
          ,
          id( nextID++ )
#endif
    {
    }
    ~Vertex()
    {
        assert( !usePtr );
    }
    BSPTreePoly::VertexUse* usePtr;
    int markVal;
#ifdef DEBUG_IDS
    int id;
    static int nextID;
#endif
};

struct BSPTreePoly::VertexUse
{
    VertexUse( Edge* edge, Vertex* vtx );
    ~VertexUse();

    void set_vertex( BSPTreePoly::Vertex*& vtx_ptr );

    BSPTreePoly::VertexUse *nextPtr, *prevPtr;
    BSPTreePoly::Vertex* vtxPtr;
    BSPTreePoly::Edge* edgePtr;
};

struct BSPTreePoly::EdgeUse
{
    EdgeUse( Edge* edge );
    EdgeUse( Edge* edge, Face* face );
    ~EdgeUse();

    BSPTreePoly::EdgeUse *prevPtr, *nextPtr;
    BSPTreePoly::Edge* edgePtr;
    BSPTreePoly::Face* facePtr;

    inline BSPTreePoly::Vertex* start() const;
    inline BSPTreePoly::Vertex* end() const;
    int sense() const;

    void insert_after( BSPTreePoly::EdgeUse* prev );
    void insert_before( BSPTreePoly::EdgeUse* next );
};

struct BSPTreePoly::Edge
{
    BSPTreePoly::VertexUse *startPtr, *endPtr;
    BSPTreePoly::EdgeUse *forwardPtr, *reversePtr;
#ifdef DEBUG_IDS
    int id;
    static int nextID;
#endif

    Edge( Vertex* vstart, Vertex* vend )
        : forwardPtr( 0 ), reversePtr( 0 )
#ifdef DEBUG_IDS
          ,
          id( nextID++ )
#endif
    {
        startPtr = new VertexUse( this, vstart );
        endPtr   = new VertexUse( this, vend );
    }

    ~Edge();

    BSPTreePoly::Vertex* start() const
    {
        return startPtr->vtxPtr;
    }
    BSPTreePoly::Vertex* end() const
    {
        return endPtr->vtxPtr;
    }

    BSPTreePoly::Face* forward() const
    {
        return forwardPtr ? forwardPtr->facePtr : 0;
    }
    BSPTreePoly::Face* reverse() const
    {
        return reversePtr ? reversePtr->facePtr : 0;
    }

    BSPTreePoly::VertexUse* use( BSPTreePoly::Vertex* vtx ) const
    {
        return ( vtx == startPtr->vtxPtr ) ? startPtr : ( vtx == endPtr->vtxPtr ) ? endPtr : 0;
    }
    BSPTreePoly::Edge* next( BSPTreePoly::Vertex* about ) const
    {
        return use( about )->nextPtr->edgePtr;
    }
    BSPTreePoly::Edge* prev( BSPTreePoly::Vertex* about ) const
    {
        return use( about )->prevPtr->edgePtr;
    }

    BSPTreePoly::EdgeUse* use( BSPTreePoly::Face* face ) const
    {
        return ( face == forwardPtr->facePtr ) ? forwardPtr : ( face == reversePtr->facePtr ) ? reversePtr : 0;
    }
    BSPTreePoly::Edge* next( BSPTreePoly::Face* about ) const
    {
        return use( about )->nextPtr->edgePtr;
    }
    BSPTreePoly::Edge* prev( BSPTreePoly::Face* about ) const
    {
        return use( about )->prevPtr->edgePtr;
    }

    BSPTreePoly::VertexUse* other( BSPTreePoly::VertexUse* vuse ) const
    {
        return vuse == startPtr ? endPtr : vuse == endPtr ? startPtr : 0;
    }
    BSPTreePoly::EdgeUse* other( BSPTreePoly::EdgeUse* vuse ) const
    {
        return vuse == forwardPtr ? reversePtr : vuse == reversePtr ? forwardPtr : 0;
    }
    BSPTreePoly::Vertex* other( BSPTreePoly::Vertex* vtx ) const
    {
        return vtx == startPtr->vtxPtr ? endPtr->vtxPtr : vtx == endPtr->vtxPtr ? startPtr->vtxPtr : 0;
    }
    BSPTreePoly::Vertex* common( BSPTreePoly::Edge* eother ) const
    {
        return start() == eother->start() || start() == eother->end() ? start()
               : end() == eother->start() || end() == eother->end()   ? end()
                                                                      : 0;
    }

    int sense( BSPTreePoly::Face* face ) const;

    void remove_from_vertex( BSPTreePoly::Vertex*& vtx_ptr );
    void remove_from_face( BSPTreePoly::Face*& face_ptr );
    void add_to_vertex( BSPTreePoly::Vertex* vtx_ptr );
};

struct BSPTreePoly::Face
{
    Face( Face* next )
        : usePtr( 0 ), nextPtr( next )
#ifdef DEBUG_IDS
          ,
          id( nextID++ )
#endif
    {
    }
    Face()
        : usePtr( 0 ), nextPtr( 0 )
#ifdef DEBUG_IDS
          ,
          id( nextID++ )
#endif
    {
    }
    ~Face();
    BSPTreePoly::EdgeUse* usePtr;
    BSPTreePoly::Face* nextPtr;
#ifdef DEBUG_IDS
    int id;
    static int nextID;
#endif
    double signed_volume() const;
};

#ifdef DEBUG_IDS
int BSPTreePoly::Vertex::nextID = 1;
int BSPTreePoly::Edge::nextID   = 1;
int BSPTreePoly::Face::nextID   = 1;
#endif
void BSPTreePoly::reset_debug_ids()
{
#ifdef DEBUG_IDS
    BSPTreePoly::Vertex::nextID = 1;
    BSPTreePoly::Edge::nextID   = 1;
    BSPTreePoly::Face::nextID   = 1;
#endif
}

// static void merge_edges( BSPTreePoly::Edge* keep_edge,
//                         BSPTreePoly::Edge* dead_edge );

static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx, BSPTreePoly::Edge* into_edge );

BSPTreePoly::VertexUse::VertexUse( BSPTreePoly::Edge* edge, BSPTreePoly::Vertex* vtx ) : vtxPtr( vtx ), edgePtr( edge )
{
    if( !vtx->usePtr )
    {
        vtx->usePtr = prevPtr = nextPtr = this;
        return;
    }

    nextPtr = vtx->usePtr;
    prevPtr = nextPtr->prevPtr;
    assert( prevPtr->nextPtr == nextPtr );
    nextPtr->prevPtr = this;
    prevPtr->nextPtr = this;
}

BSPTreePoly::VertexUse::~VertexUse()
{
    if( nextPtr == this )
    {
        assert( prevPtr == this );
        assert( vtxPtr->usePtr == this );
        vtxPtr->usePtr = 0;
        delete vtxPtr;
    }
    else if( vtxPtr->usePtr == this )
        vtxPtr->usePtr = nextPtr;

    nextPtr->prevPtr = prevPtr;
    prevPtr->nextPtr = nextPtr;
    nextPtr = prevPtr = 0;
}

void BSPTreePoly::VertexUse::set_vertex( BSPTreePoly::Vertex*& vtx )
{
    if( vtxPtr )
    {
        if( nextPtr == prevPtr )
        {
            assert( nextPtr == this );
            vtxPtr->usePtr = 0;
            delete vtx;
            vtx = 0;
        }
        else
        {
            nextPtr->prevPtr = prevPtr;
            prevPtr->nextPtr = nextPtr;
            if( vtxPtr->usePtr == this ) vtxPtr->usePtr = nextPtr;
        }
    }

    if( vtx )
    {
        vtxPtr                  = vtx;
        nextPtr                 = vtxPtr->usePtr->nextPtr;
        prevPtr                 = vtxPtr->usePtr;
        nextPtr->prevPtr        = this;
        vtxPtr->usePtr->nextPtr = this;
    }
}

BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge ) : prevPtr( 0 ), nextPtr( 0 ), edgePtr( edge ), facePtr( 0 ) {}

BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge, BSPTreePoly::Face* face ) : edgePtr( edge ), facePtr( face )
{
    assert( !face->usePtr );
    face->usePtr = prevPtr = nextPtr = this;

    if( !face->usePtr )
    {
        face->usePtr = prevPtr = nextPtr = this;
        return;
    }

    nextPtr = face->usePtr;
    prevPtr = nextPtr->prevPtr;
    assert( prevPtr->nextPtr == nextPtr );
    nextPtr->prevPtr = this;
    prevPtr->nextPtr = this;
}

void BSPTreePoly::EdgeUse::insert_after( BSPTreePoly::EdgeUse* prev )
{
    // shouldn't already be in a face
    assert( !facePtr );
    // adjacent edges should share vertices
    assert( start() == prev->end() );

    facePtr          = prev->facePtr;
    nextPtr          = prev->nextPtr;
    prevPtr          = prev;
    nextPtr->prevPtr = this;
    prevPtr->nextPtr = this;
}

void BSPTreePoly::EdgeUse::insert_before( BSPTreePoly::EdgeUse* next )
{
    // shouldn't already be in a face
    assert( !facePtr );
    // adjacent edges should share vertices
    assert( end() == next->start() );

    facePtr          = next->facePtr;
    prevPtr          = next->prevPtr;
    nextPtr          = next;
    nextPtr->prevPtr = this;
    prevPtr->nextPtr = this;
}

BSPTreePoly::EdgeUse::~EdgeUse()
{
    if( facePtr->usePtr == this ) facePtr->usePtr = ( nextPtr == this ) ? 0 : nextPtr;

    if( edgePtr->forwardPtr == this ) edgePtr->forwardPtr = 0;
    if( edgePtr->reversePtr == this ) edgePtr->reversePtr = 0;

    if( !edgePtr->forwardPtr && !edgePtr->reversePtr ) delete edgePtr;

    nextPtr->prevPtr = prevPtr;
    prevPtr->nextPtr = nextPtr;
    nextPtr = prevPtr = 0;
}

int BSPTreePoly::EdgeUse::sense() const
{
    if( edgePtr->forwardPtr == this )
        return 1;
    else if( edgePtr->reversePtr == this )
        return -1;
    else
        return 0;
}

BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::start() const
{
    if( edgePtr->forwardPtr == this )
        return edgePtr->start();
    else if( edgePtr->reversePtr == this )
        return edgePtr->end();
    else
        return 0;
}

BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::end() const
{
    if( edgePtr->forwardPtr == this )
        return edgePtr->end();
    else if( edgePtr->reversePtr == this )
        return edgePtr->start();
    else
        return 0;
}

BSPTreePoly::Edge::~Edge()
{
    delete startPtr;
    delete endPtr;
    delete forwardPtr;
    delete reversePtr;
}

int BSPTreePoly::Edge::sense( BSPTreePoly::Face* face ) const
{
    if( forwardPtr && forwardPtr->facePtr == face )
        return 1;
    else if( reversePtr && reversePtr->facePtr == face )
        return -1;
    else
        return 0;
}

BSPTreePoly::Face::~Face()
{
    BSPTreePoly::EdgeUse* nextEdgeUsePtr = usePtr;
    while( nextEdgeUsePtr )
    {
        delete nextEdgeUsePtr;  // This is tricky: ~EdgeUse() might change the value of usePtr
        if( usePtr && usePtr != nextEdgeUsePtr )
            nextEdgeUsePtr = usePtr;
        else
            nextEdgeUsePtr = 0;
    }
    usePtr = 0;
}

void BSPTreePoly::clear()
{
    while( faceList )
    {
        Face* face = faceList;
        faceList   = faceList->nextPtr;
        delete face;
    }
}

ErrorCode BSPTreePoly::set( const CartVect hex_corners[8] )
{
    clear();

    Vertex* vertices[8];
    for( int i = 0; i < 8; ++i )
        vertices[i] = new Vertex( hex_corners[i] );

    Edge* edges[12];
#ifdef DEBUG_IDS
    int start_id = Edge::nextID;
#endif
    for( int i = 0; i < 4; ++i )
    {
        int j        = ( i + 1 ) % 4;
        edges[i]     = new Edge( vertices[i], vertices[j] );
        edges[i + 4] = new Edge( vertices[i], vertices[i + 4] );
        edges[i + 8] = new Edge( vertices[i + 4], vertices[j + 4] );
    }
#ifdef DEBUG_IDS
    for( int i = 0; i < 12; ++i )
        edges[i]->id = start_id++;
#endif

    static const int face_conn[6][4] = { { 0, 5, -8, -4 },  { 1, 6, -9, -5 },    { 2, 7, -10, -6 },
                                         { 3, 4, -11, -7 }, { -3, -2, -1, -12 }, { 8, 9, 10, 11 } };
    for( int i = 0; i < 6; ++i )
    {
        faceList      = new Face( faceList );
        EdgeUse* prev = 0;
        for( int j = 0; j < 4; ++j )
        {
            int e = face_conn[i][j];
            if( e < 0 )
            {
                e = ( -e ) % 12;
                assert( !edges[e]->reversePtr );
                if( !prev )
                {
                    edges[e]->reversePtr = new EdgeUse( edges[e], faceList );
                }
                else
                {
                    edges[e]->reversePtr = new EdgeUse( edges[e] );
                    edges[e]->reversePtr->insert_after( prev );
                }
                prev = edges[e]->reversePtr;
            }
            else
            {
                assert( !edges[e]->forwardPtr );
                if( !prev )
                {
                    edges[e]->forwardPtr = new EdgeUse( edges[e], faceList );
                }
                else
                {
                    edges[e]->forwardPtr = new EdgeUse( edges[e] );
                    edges[e]->forwardPtr->insert_after( prev );
                }
                prev = edges[e]->forwardPtr;
            }
        }
    }

    return MB_SUCCESS;
}

void BSPTreePoly::get_faces( std::vector< const Face* >& face_list ) const
{
    face_list.clear();
    for( Face* face = faceList; face; face = face->nextPtr )
        face_list.push_back( face );
}

void BSPTreePoly::get_vertices( const Face* face, std::vector< CartVect >& vertices ) const
{
    vertices.clear();
    if( !face || !face->usePtr ) return;

    EdgeUse* coedge = face->usePtr;
    do
    {
        vertices.push_back( *coedge->end() );
        coedge = coedge->nextPtr;
    } while( coedge != face->usePtr );
}

double BSPTreePoly::Face::signed_volume() const
{
    CartVect sum( 0.0 );
    const CartVect* base = usePtr->start();
    CartVect d1          = ( *usePtr->end() - *base );
    for( EdgeUse* coedge = usePtr->nextPtr; coedge != usePtr; coedge = coedge->nextPtr )
    {
        CartVect d2 = ( *coedge->end() - *base );
        sum += d1 * d2;
        d1 = d2;
    }
    return ( 1.0 / 6.0 ) * ( sum % *base );
}

double BSPTreePoly::volume() const
{
    double result = 0;
    for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
        result += ptr->signed_volume();
    return result;
}

void BSPTreePoly::set_vertex_marks( int value )
{
    for( Face* face = faceList; face; face = face->nextPtr )
    {
        EdgeUse* edge = face->usePtr;
        do
        {
            edge->edgePtr->start()->markVal = value;
            edge->edgePtr->end()->markVal   = value;
            edge                            = edge->nextPtr;
        } while( edge && edge != face->usePtr );
    }
}
/*
static void merge_edges( BSPTreePoly::Edge* keep_edge,
                         BSPTreePoly::Edge* dead_edge )
{
  // edges must share a vertex
  BSPTreePoly::Vertex* dead_vtx = keep_edge->common(dead_edge);
  assert(dead_vtx);
   // vertex may have only two adjacent edges
  BSPTreePoly::VertexUse* dead_vtxuse = dead_edge->use(dead_vtx);
  assert(dead_vtxuse);
  BSPTreePoly::VertexUse* keep_vtxuse = dead_vtxuse->nextPtr;
  assert(keep_vtxuse);
  assert(keep_vtxuse->edgePtr == keep_edge);
  assert(keep_vtxuse->nextPtr == dead_vtxuse);
  assert(keep_vtxuse->prevPtr == dead_vtxuse);
  assert(dead_vtxuse->prevPtr == keep_vtxuse);

  // kept edge now ends with the kept vertex on the dead edge
  keep_vtxuse->set_vertex( dead_edge->other(dead_vtx) );

  // destructors should take care of everything else
  // (including removing dead edge from face loops)
  delete dead_edge;
}
*/
static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx, BSPTreePoly::Edge* into_edge )
{
    // split edge, creating new edge
    BSPTreePoly::Edge* new_edge = new BSPTreePoly::Edge( new_vtx, into_edge->end() );
    into_edge->endPtr->set_vertex( new_vtx );  // This call might delete new_vtx

    // update coedge loops in faces
    if( into_edge->forwardPtr )
    {
        new_edge->forwardPtr = new BSPTreePoly::EdgeUse( new_edge );
        new_edge->forwardPtr->insert_after( into_edge->forwardPtr );
    }
    if( into_edge->reversePtr )
    {
        new_edge->reversePtr = new BSPTreePoly::EdgeUse( new_edge );
        new_edge->reversePtr->insert_before( into_edge->reversePtr );
    }

    return new_edge;
}

static BSPTreePoly::Face* split_face( BSPTreePoly::EdgeUse* start, BSPTreePoly::EdgeUse* end )
{
    BSPTreePoly::Face* face = start->facePtr;
    assert( face == end->facePtr );
    BSPTreePoly::Face* new_face      = new BSPTreePoly::Face;
    BSPTreePoly::EdgeUse* keep_start = start->prevPtr;
    BSPTreePoly::EdgeUse* keep_end   = end->nextPtr;
    for( BSPTreePoly::EdgeUse* ptr = start; ptr != keep_end; ptr = ptr->nextPtr )
    {
        if( face->usePtr == ptr ) face->usePtr = keep_start;
        ptr->facePtr = new_face;
    }
    new_face->usePtr        = start;
    BSPTreePoly::Edge* edge = new BSPTreePoly::Edge( start->start(), end->end() );
    edge->forwardPtr        = new BSPTreePoly::EdgeUse( edge );
    edge->reversePtr        = new BSPTreePoly::EdgeUse( edge );

    edge->forwardPtr->facePtr = face;
    edge->forwardPtr->prevPtr = keep_start;
    keep_start->nextPtr       = edge->forwardPtr;
    edge->forwardPtr->nextPtr = keep_end;
    keep_end->prevPtr         = edge->forwardPtr;

    edge->reversePtr->facePtr = new_face;
    edge->reversePtr->nextPtr = start;
    start->prevPtr            = edge->reversePtr;
    edge->reversePtr->prevPtr = end;
    end->nextPtr              = edge->reversePtr;

    return new_face;
}

bool BSPTreePoly::cut_polyhedron( const CartVect& plane_normal, double plane_coeff )
{
    const double EPSILON = 1e-6;  // points this close are considered coincident

    // scale epsilon rather than normalizing normal vector
    const double epsilon = EPSILON * ( plane_normal % plane_normal );

    // Classify all points above/below plane and destroy any faces
    // that have no vertices below the plane.
    const int UNKNOWN = 0;
    const int ABOVE   = 1;
    const int ON      = 2;
    const int BELOW   = 3;
    int num_above     = 0;
    set_vertex_marks( UNKNOWN );

    // Classify all points above/below plane and
    // split any edge that intersect the plane.
    for( Face* face = faceList; face; face = face->nextPtr )
    {
        EdgeUse* edge = face->usePtr;

        do
        {
            Vertex* start = edge->edgePtr->start();
            Vertex* end   = edge->edgePtr->end();

            if( !start->markVal )
            {
                double d = plane_normal % *start + plane_coeff;
                if( d * d <= epsilon )
                    start->markVal = ON;
                else if( d < 0.0 )
                    start->markVal = BELOW;
                else
                {
                    start->markVal = ABOVE;
                    ++num_above;
                }
            }

            if( !end->markVal )
            {
                double d = plane_normal % *end + plane_coeff;
                if( d * d <= epsilon )
                    end->markVal = ON;
                else if( d < 0.0 )
                    end->markVal = BELOW;
                else
                {
                    end->markVal = ABOVE;
                    ++num_above;
                }
            }

            if( ( end->markVal == ABOVE && start->markVal == BELOW ) ||
                ( end->markVal == BELOW && start->markVal == ABOVE ) )
            {
                CartVect dir     = *end - *start;
                double t         = -( plane_normal % *start + plane_coeff ) / ( dir % plane_normal );
                Vertex* new_vtx  = new Vertex( *start + t * dir );
                new_vtx->markVal = ON;
                split_edge( new_vtx, edge->edgePtr );  // This call might delete new_vtx
                end = new_vtx;
            }

            edge = edge->nextPtr;
        } while( edge && edge != face->usePtr );
    }

    if( !num_above ) return false;

    // Split faces
    for( Face* face = faceList; face; face = face->nextPtr )
    {
        EdgeUse* edge = face->usePtr;

        EdgeUse *split_start = 0, *split_end = 0, *other_split = 0;
        do
        {
            if( edge->end()->markVal == ON && edge->start()->markVal != ON )
            {
                if( !split_start )
                    split_start = edge->nextPtr;
                else if( !split_end )
                    split_end = edge;
                else
                    other_split = edge;
            }

            edge = edge->nextPtr;
        } while( edge && edge != face->usePtr );

        // If two vertices are on plane (but not every vertex)
        // then split the face
        if( split_end && !other_split )
        {
            assert( split_start );
            Face* new_face    = split_face( split_start, split_end );
            new_face->nextPtr = faceList;
            faceList          = new_face;
        }
    }

    // Destroy all faces that are above the plane
    Face** lptr = &faceList;
    while( *lptr )
    {
        EdgeUse* edge   = ( *lptr )->usePtr;
        bool some_above = false;
        do
        {
            if( edge->start()->markVal == ABOVE )
            {
                some_above = true;
                break;
            }
            edge = edge->nextPtr;
        } while( edge && edge != ( *lptr )->usePtr );

        if( some_above )
        {
            Face* dead = *lptr;
            *lptr      = ( *lptr )->nextPtr;
            delete dead;
        }
        else
        {
            lptr = &( ( *lptr )->nextPtr );
        }
    }

    // Construct a new face in the cut plane

    // First find an edge to start at
    Edge* edge_ptr = 0;
    for( Face* face = faceList; face && !edge_ptr; face = face->nextPtr )
    {
        EdgeUse* co_edge = face->usePtr;
        do
        {
            if( 0 == co_edge->edgePtr->other( co_edge ) )
            {
                edge_ptr = co_edge->edgePtr;
                break;
            }
            co_edge = co_edge->nextPtr;
        } while( co_edge && co_edge != face->usePtr );
    }
    if( !edge_ptr ) return false;

    // Constuct new face and first CoEdge
    faceList = new Face( faceList );
    Vertex *next_vtx, *start_vtx;
    EdgeUse* prev_coedge;
    if( edge_ptr->forwardPtr )
    {
        next_vtx  = edge_ptr->start();
        start_vtx = edge_ptr->end();
        assert( !edge_ptr->reversePtr );
        prev_coedge = edge_ptr->reversePtr = new EdgeUse( edge_ptr, faceList );
    }
    else
    {
        next_vtx    = edge_ptr->end();
        start_vtx   = edge_ptr->start();
        prev_coedge = edge_ptr->forwardPtr = new EdgeUse( edge_ptr, faceList );
    }

    // Construct coedges until loop is closed
    while( next_vtx != start_vtx )
    {
        // find next edge adjacent to vertex with only one adjacent face
        VertexUse* this_use = edge_ptr->use( next_vtx );
        VertexUse* use      = this_use->nextPtr;
        while( use != this_use )
        {
            if( use->edgePtr->forwardPtr == 0 )
            {
                edge_ptr = use->edgePtr;
                assert( edge_ptr->start() == next_vtx );
                next_vtx             = edge_ptr->end();
                edge_ptr->forwardPtr = new EdgeUse( edge_ptr );
                edge_ptr->forwardPtr->insert_after( prev_coedge );
                prev_coedge = edge_ptr->forwardPtr;
                break;
            }
            else if( use->edgePtr->reversePtr == 0 )
            {
                edge_ptr = use->edgePtr;
                assert( edge_ptr->end() == next_vtx );
                next_vtx             = edge_ptr->start();
                edge_ptr->reversePtr = new EdgeUse( edge_ptr );
                edge_ptr->reversePtr->insert_after( prev_coedge );
                prev_coedge = edge_ptr->reversePtr;
                break;
            }

            use = use->nextPtr;
            assert( use != this_use );  // failed to close loop!
        }
    }

    return true;
}

bool BSPTreePoly::is_valid() const
{
    std::set< Face* > list_faces;

    int i = 0;
    for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
    {
        if( ++i > 10000 ) return false;
        if( !list_faces.insert( ptr ).second ) return false;
    }

    std::set< Vertex* > vertices;
    for( Face* face = faceList; face; face = face->nextPtr )
    {
        i               = 0;
        EdgeUse* coedge = face->usePtr;
        do
        {
            if( ++i > 10000 ) return false;

            if( coedge->facePtr != face ) return false;

            Edge* edge = coedge->edgePtr;
            if( !edge->startPtr || !edge->endPtr ) return false;

            vertices.insert( edge->start() );
            vertices.insert( edge->end() );

            EdgeUse* other;
            if( edge->forwardPtr == coedge )
                other = edge->reversePtr;
            else if( edge->reversePtr != coedge )
                return false;
            else
                other = edge->forwardPtr;
            if( !other ) return false;
            if( list_faces.find( other->facePtr ) == list_faces.end() ) return false;

            EdgeUse* next = coedge->nextPtr;
            if( next->prevPtr != coedge ) return false;
            if( coedge->end() != next->start() ) return false;

            coedge = next;
        } while( coedge != face->usePtr );
    }

    for( std::set< Vertex* >::iterator j = vertices.begin(); j != vertices.end(); ++j )
    {
        Vertex* vtx = *j;

        i              = 0;
        VertexUse* use = vtx->usePtr;
        do
        {
            if( ++i > 10000 ) return false;

            if( use->vtxPtr != vtx ) return false;

            Edge* edge = use->edgePtr;
            if( !edge ) return false;
            if( edge->startPtr != use && edge->endPtr != use ) return false;

            VertexUse* next = use->nextPtr;
            if( next->prevPtr != use ) return false;

            use = next;
        } while( use != vtx->usePtr );
    }

    return true;
}

bool BSPTreePoly::is_point_contained( const CartVect& point ) const
{
    if( !faceList )  // empty (zero-dimension) polyhedron
        return false;

    const double EPSILON = 1e-6;
    // Test that point is below the plane of each face
    // NOTE: This will NOT work for polyhedra w/ concavities
    for( Face* face = faceList; face; face = face->nextPtr )
    {
        Vertex *pt1, *pt2, *pt3;
        pt1 = face->usePtr->start();
        pt2 = face->usePtr->end();
        pt3 = face->usePtr->nextPtr->end();

        if( pt3 == pt1 )  // degenerate
            continue;

        CartVect norm = ( *pt3 - *pt2 ) * ( *pt1 - *pt2 );
        double coeff  = -( norm % *pt2 );
        if( ( norm % point + coeff ) > EPSILON )  // if above plane, with some -epsilon
            return false;
    }

    return true;
}

}  // namespace moab