MOAB: Mesh Oriented datABase  (version 5.4.1)
BSPTreePoly.cpp
Go to the documentation of this file.
00001 #include "moab/CartVect.hpp"
00002 #include "moab/BSPTreePoly.hpp"
00003 #include <cassert>
00004 #include <cstdlib>
00005 #include <set>
00006 
00007 #undef DEBUG_IDS
00008 
00009 namespace moab
00010 {
00011 
00012 struct BSPTreePoly::Vertex : public CartVect
00013 {
00014     Vertex( const CartVect& v )
00015         : CartVect( v ), usePtr( 0 ), markVal( 0 )
00016 #ifdef DEBUG_IDS
00017           ,
00018           id( nextID++ )
00019 #endif
00020     {
00021     }
00022     ~Vertex()
00023     {
00024         assert( !usePtr );
00025     }
00026     BSPTreePoly::VertexUse* usePtr;
00027     int markVal;
00028 #ifdef DEBUG_IDS
00029     int id;
00030     static int nextID;
00031 #endif
00032 };
00033 
00034 struct BSPTreePoly::VertexUse
00035 {
00036     VertexUse( Edge* edge, Vertex* vtx );
00037     ~VertexUse();
00038 
00039     void set_vertex( BSPTreePoly::Vertex*& vtx_ptr );
00040 
00041     BSPTreePoly::VertexUse *nextPtr, *prevPtr;
00042     BSPTreePoly::Vertex* vtxPtr;
00043     BSPTreePoly::Edge* edgePtr;
00044 };
00045 
00046 struct BSPTreePoly::EdgeUse
00047 {
00048     EdgeUse( Edge* edge );
00049     EdgeUse( Edge* edge, Face* face );
00050     ~EdgeUse();
00051 
00052     BSPTreePoly::EdgeUse *prevPtr, *nextPtr;
00053     BSPTreePoly::Edge* edgePtr;
00054     BSPTreePoly::Face* facePtr;
00055 
00056     inline BSPTreePoly::Vertex* start() const;
00057     inline BSPTreePoly::Vertex* end() const;
00058     int sense() const;
00059 
00060     void insert_after( BSPTreePoly::EdgeUse* prev );
00061     void insert_before( BSPTreePoly::EdgeUse* next );
00062 };
00063 
00064 struct BSPTreePoly::Edge
00065 {
00066     BSPTreePoly::VertexUse *startPtr, *endPtr;
00067     BSPTreePoly::EdgeUse *forwardPtr, *reversePtr;
00068 #ifdef DEBUG_IDS
00069     int id;
00070     static int nextID;
00071 #endif
00072 
00073     Edge( Vertex* vstart, Vertex* vend )
00074         : forwardPtr( 0 ), reversePtr( 0 )
00075 #ifdef DEBUG_IDS
00076           ,
00077           id( nextID++ )
00078 #endif
00079     {
00080         startPtr = new VertexUse( this, vstart );
00081         endPtr   = new VertexUse( this, vend );
00082     }
00083 
00084     ~Edge();
00085 
00086     BSPTreePoly::Vertex* start() const
00087     {
00088         return startPtr->vtxPtr;
00089     }
00090     BSPTreePoly::Vertex* end() const
00091     {
00092         return endPtr->vtxPtr;
00093     }
00094 
00095     BSPTreePoly::Face* forward() const
00096     {
00097         return forwardPtr ? forwardPtr->facePtr : 0;
00098     }
00099     BSPTreePoly::Face* reverse() const
00100     {
00101         return reversePtr ? reversePtr->facePtr : 0;
00102     }
00103 
00104     BSPTreePoly::VertexUse* use( BSPTreePoly::Vertex* vtx ) const
00105     {
00106         return ( vtx == startPtr->vtxPtr ) ? startPtr : ( vtx == endPtr->vtxPtr ) ? endPtr : 0;
00107     }
00108     BSPTreePoly::Edge* next( BSPTreePoly::Vertex* about ) const
00109     {
00110         return use( about )->nextPtr->edgePtr;
00111     }
00112     BSPTreePoly::Edge* prev( BSPTreePoly::Vertex* about ) const
00113     {
00114         return use( about )->prevPtr->edgePtr;
00115     }
00116 
00117     BSPTreePoly::EdgeUse* use( BSPTreePoly::Face* face ) const
00118     {
00119         return ( face == forwardPtr->facePtr ) ? forwardPtr : ( face == reversePtr->facePtr ) ? reversePtr : 0;
00120     }
00121     BSPTreePoly::Edge* next( BSPTreePoly::Face* about ) const
00122     {
00123         return use( about )->nextPtr->edgePtr;
00124     }
00125     BSPTreePoly::Edge* prev( BSPTreePoly::Face* about ) const
00126     {
00127         return use( about )->prevPtr->edgePtr;
00128     }
00129 
00130     BSPTreePoly::VertexUse* other( BSPTreePoly::VertexUse* vuse ) const
00131     {
00132         return vuse == startPtr ? endPtr : vuse == endPtr ? startPtr : 0;
00133     }
00134     BSPTreePoly::EdgeUse* other( BSPTreePoly::EdgeUse* vuse ) const
00135     {
00136         return vuse == forwardPtr ? reversePtr : vuse == reversePtr ? forwardPtr : 0;
00137     }
00138     BSPTreePoly::Vertex* other( BSPTreePoly::Vertex* vtx ) const
00139     {
00140         return vtx == startPtr->vtxPtr ? endPtr->vtxPtr : vtx == endPtr->vtxPtr ? startPtr->vtxPtr : 0;
00141     }
00142     BSPTreePoly::Vertex* common( BSPTreePoly::Edge* eother ) const
00143     {
00144         return start() == eother->start() || start() == eother->end() ? start()
00145                : end() == eother->start() || end() == eother->end()   ? end()
00146                                                                       : 0;
00147     }
00148 
00149     int sense( BSPTreePoly::Face* face ) const;
00150 
00151     void remove_from_vertex( BSPTreePoly::Vertex*& vtx_ptr );
00152     void remove_from_face( BSPTreePoly::Face*& face_ptr );
00153     void add_to_vertex( BSPTreePoly::Vertex* vtx_ptr );
00154 };
00155 
00156 struct BSPTreePoly::Face
00157 {
00158     Face( Face* next )
00159         : usePtr( 0 ), nextPtr( next )
00160 #ifdef DEBUG_IDS
00161           ,
00162           id( nextID++ )
00163 #endif
00164     {
00165     }
00166     Face()
00167         : usePtr( 0 ), nextPtr( 0 )
00168 #ifdef DEBUG_IDS
00169           ,
00170           id( nextID++ )
00171 #endif
00172     {
00173     }
00174     ~Face();
00175     BSPTreePoly::EdgeUse* usePtr;
00176     BSPTreePoly::Face* nextPtr;
00177 #ifdef DEBUG_IDS
00178     int id;
00179     static int nextID;
00180 #endif
00181     double signed_volume() const;
00182 };
00183 
00184 #ifdef DEBUG_IDS
00185 int BSPTreePoly::Vertex::nextID = 1;
00186 int BSPTreePoly::Edge::nextID   = 1;
00187 int BSPTreePoly::Face::nextID   = 1;
00188 #endif
00189 void BSPTreePoly::reset_debug_ids()
00190 {
00191 #ifdef DEBUG_IDS
00192     BSPTreePoly::Vertex::nextID = 1;
00193     BSPTreePoly::Edge::nextID   = 1;
00194     BSPTreePoly::Face::nextID   = 1;
00195 #endif
00196 }
00197 
00198 // static void merge_edges( BSPTreePoly::Edge* keep_edge,
00199 //                         BSPTreePoly::Edge* dead_edge );
00200 
00201 static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx, BSPTreePoly::Edge* into_edge );
00202 
00203 BSPTreePoly::VertexUse::VertexUse( BSPTreePoly::Edge* edge, BSPTreePoly::Vertex* vtx ) : vtxPtr( vtx ), edgePtr( edge )
00204 {
00205     if( !vtx->usePtr )
00206     {
00207         vtx->usePtr = prevPtr = nextPtr = this;
00208         return;
00209     }
00210 
00211     nextPtr = vtx->usePtr;
00212     prevPtr = nextPtr->prevPtr;
00213     assert( prevPtr->nextPtr == nextPtr );
00214     nextPtr->prevPtr = this;
00215     prevPtr->nextPtr = this;
00216 }
00217 
00218 BSPTreePoly::VertexUse::~VertexUse()
00219 {
00220     if( nextPtr == this )
00221     {
00222         assert( prevPtr == this );
00223         assert( vtxPtr->usePtr == this );
00224         vtxPtr->usePtr = 0;
00225         delete vtxPtr;
00226     }
00227     else if( vtxPtr->usePtr == this )
00228         vtxPtr->usePtr = nextPtr;
00229 
00230     nextPtr->prevPtr = prevPtr;
00231     prevPtr->nextPtr = nextPtr;
00232     nextPtr = prevPtr = 0;
00233 }
00234 
00235 void BSPTreePoly::VertexUse::set_vertex( BSPTreePoly::Vertex*& vtx )
00236 {
00237     if( vtxPtr )
00238     {
00239         if( nextPtr == prevPtr )
00240         {
00241             assert( nextPtr == this );
00242             vtxPtr->usePtr = 0;
00243             delete vtx;
00244             vtx = 0;
00245         }
00246         else
00247         {
00248             nextPtr->prevPtr = prevPtr;
00249             prevPtr->nextPtr = nextPtr;
00250             if( vtxPtr->usePtr == this ) vtxPtr->usePtr = nextPtr;
00251         }
00252     }
00253 
00254     if( vtx )
00255     {
00256         vtxPtr                  = vtx;
00257         nextPtr                 = vtxPtr->usePtr->nextPtr;
00258         prevPtr                 = vtxPtr->usePtr;
00259         nextPtr->prevPtr        = this;
00260         vtxPtr->usePtr->nextPtr = this;
00261     }
00262 }
00263 
00264 BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge ) : prevPtr( 0 ), nextPtr( 0 ), edgePtr( edge ), facePtr( 0 ) {}
00265 
00266 BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge, BSPTreePoly::Face* face ) : edgePtr( edge ), facePtr( face )
00267 {
00268     assert( !face->usePtr );
00269     face->usePtr = prevPtr = nextPtr = this;
00270 
00271     if( !face->usePtr )
00272     {
00273         face->usePtr = prevPtr = nextPtr = this;
00274         return;
00275     }
00276 
00277     nextPtr = face->usePtr;
00278     prevPtr = nextPtr->prevPtr;
00279     assert( prevPtr->nextPtr == nextPtr );
00280     nextPtr->prevPtr = this;
00281     prevPtr->nextPtr = this;
00282 }
00283 
00284 void BSPTreePoly::EdgeUse::insert_after( BSPTreePoly::EdgeUse* prev )
00285 {
00286     // shouldn't already be in a face
00287     assert( !facePtr );
00288     // adjacent edges should share vertices
00289     assert( start() == prev->end() );
00290 
00291     facePtr          = prev->facePtr;
00292     nextPtr          = prev->nextPtr;
00293     prevPtr          = prev;
00294     nextPtr->prevPtr = this;
00295     prevPtr->nextPtr = this;
00296 }
00297 
00298 void BSPTreePoly::EdgeUse::insert_before( BSPTreePoly::EdgeUse* next )
00299 {
00300     // shouldn't already be in a face
00301     assert( !facePtr );
00302     // adjacent edges should share vertices
00303     assert( end() == next->start() );
00304 
00305     facePtr          = next->facePtr;
00306     prevPtr          = next->prevPtr;
00307     nextPtr          = next;
00308     nextPtr->prevPtr = this;
00309     prevPtr->nextPtr = this;
00310 }
00311 
00312 BSPTreePoly::EdgeUse::~EdgeUse()
00313 {
00314     if( facePtr->usePtr == this ) facePtr->usePtr = ( nextPtr == this ) ? 0 : nextPtr;
00315 
00316     if( edgePtr->forwardPtr == this ) edgePtr->forwardPtr = 0;
00317     if( edgePtr->reversePtr == this ) edgePtr->reversePtr = 0;
00318 
00319     if( !edgePtr->forwardPtr && !edgePtr->reversePtr ) delete edgePtr;
00320 
00321     nextPtr->prevPtr = prevPtr;
00322     prevPtr->nextPtr = nextPtr;
00323     nextPtr = prevPtr = 0;
00324 }
00325 
00326 int BSPTreePoly::EdgeUse::sense() const
00327 {
00328     if( edgePtr->forwardPtr == this )
00329         return 1;
00330     else if( edgePtr->reversePtr == this )
00331         return -1;
00332     else
00333         return 0;
00334 }
00335 
00336 BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::start() const
00337 {
00338     if( edgePtr->forwardPtr == this )
00339         return edgePtr->start();
00340     else if( edgePtr->reversePtr == this )
00341         return edgePtr->end();
00342     else
00343         return 0;
00344 }
00345 
00346 BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::end() const
00347 {
00348     if( edgePtr->forwardPtr == this )
00349         return edgePtr->end();
00350     else if( edgePtr->reversePtr == this )
00351         return edgePtr->start();
00352     else
00353         return 0;
00354 }
00355 
00356 BSPTreePoly::Edge::~Edge()
00357 {
00358     delete startPtr;
00359     delete endPtr;
00360     delete forwardPtr;
00361     delete reversePtr;
00362 }
00363 
00364 int BSPTreePoly::Edge::sense( BSPTreePoly::Face* face ) const
00365 {
00366     if( forwardPtr && forwardPtr->facePtr == face )
00367         return 1;
00368     else if( reversePtr && reversePtr->facePtr == face )
00369         return -1;
00370     else
00371         return 0;
00372 }
00373 
00374 BSPTreePoly::Face::~Face()
00375 {
00376     BSPTreePoly::EdgeUse* nextEdgeUsePtr = usePtr;
00377     while( nextEdgeUsePtr )
00378     {
00379         delete nextEdgeUsePtr;  // This is tricky: ~EdgeUse() might change the value of usePtr
00380         if( usePtr && usePtr != nextEdgeUsePtr )
00381             nextEdgeUsePtr = usePtr;
00382         else
00383             nextEdgeUsePtr = 0;
00384     }
00385     usePtr = 0;
00386 }
00387 
00388 void BSPTreePoly::clear()
00389 {
00390     while( faceList )
00391     {
00392         Face* face = faceList;
00393         faceList   = faceList->nextPtr;
00394         delete face;
00395     }
00396 }
00397 
00398 ErrorCode BSPTreePoly::set( const CartVect hex_corners[8] )
00399 {
00400     clear();
00401 
00402     Vertex* vertices[8];
00403     for( int i = 0; i < 8; ++i )
00404         vertices[i] = new Vertex( hex_corners[i] );
00405 
00406     Edge* edges[12];
00407 #ifdef DEBUG_IDS
00408     int start_id = Edge::nextID;
00409 #endif
00410     for( int i = 0; i < 4; ++i )
00411     {
00412         int j        = ( i + 1 ) % 4;
00413         edges[i]     = new Edge( vertices[i], vertices[j] );
00414         edges[i + 4] = new Edge( vertices[i], vertices[i + 4] );
00415         edges[i + 8] = new Edge( vertices[i + 4], vertices[j + 4] );
00416     }
00417 #ifdef DEBUG_IDS
00418     for( int i = 0; i < 12; ++i )
00419         edges[i]->id = start_id++;
00420 #endif
00421 
00422     static const int face_conn[6][4] = { { 0, 5, -8, -4 },  { 1, 6, -9, -5 },    { 2, 7, -10, -6 },
00423                                          { 3, 4, -11, -7 }, { -3, -2, -1, -12 }, { 8, 9, 10, 11 } };
00424     for( int i = 0; i < 6; ++i )
00425     {
00426         faceList      = new Face( faceList );
00427         EdgeUse* prev = 0;
00428         for( int j = 0; j < 4; ++j )
00429         {
00430             int e = face_conn[i][j];
00431             if( e < 0 )
00432             {
00433                 e = ( -e ) % 12;
00434                 assert( !edges[e]->reversePtr );
00435                 if( !prev )
00436                 {
00437                     edges[e]->reversePtr = new EdgeUse( edges[e], faceList );
00438                 }
00439                 else
00440                 {
00441                     edges[e]->reversePtr = new EdgeUse( edges[e] );
00442                     edges[e]->reversePtr->insert_after( prev );
00443                 }
00444                 prev = edges[e]->reversePtr;
00445             }
00446             else
00447             {
00448                 assert( !edges[e]->forwardPtr );
00449                 if( !prev )
00450                 {
00451                     edges[e]->forwardPtr = new EdgeUse( edges[e], faceList );
00452                 }
00453                 else
00454                 {
00455                     edges[e]->forwardPtr = new EdgeUse( edges[e] );
00456                     edges[e]->forwardPtr->insert_after( prev );
00457                 }
00458                 prev = edges[e]->forwardPtr;
00459             }
00460         }
00461     }
00462 
00463     return MB_SUCCESS;
00464 }
00465 
00466 void BSPTreePoly::get_faces( std::vector< const Face* >& face_list ) const
00467 {
00468     face_list.clear();
00469     for( Face* face = faceList; face; face = face->nextPtr )
00470         face_list.push_back( face );
00471 }
00472 
00473 void BSPTreePoly::get_vertices( const Face* face, std::vector< CartVect >& vertices ) const
00474 {
00475     vertices.clear();
00476     if( !face || !face->usePtr ) return;
00477 
00478     EdgeUse* coedge = face->usePtr;
00479     do
00480     {
00481         vertices.push_back( *coedge->end() );
00482         coedge = coedge->nextPtr;
00483     } while( coedge != face->usePtr );
00484 }
00485 
00486 double BSPTreePoly::Face::signed_volume() const
00487 {
00488     CartVect sum( 0.0 );
00489     const CartVect* base = usePtr->start();
00490     CartVect d1          = ( *usePtr->end() - *base );
00491     for( EdgeUse* coedge = usePtr->nextPtr; coedge != usePtr; coedge = coedge->nextPtr )
00492     {
00493         CartVect d2 = ( *coedge->end() - *base );
00494         sum += d1 * d2;
00495         d1 = d2;
00496     }
00497     return ( 1.0 / 6.0 ) * ( sum % *base );
00498 }
00499 
00500 double BSPTreePoly::volume() const
00501 {
00502     double result = 0;
00503     for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
00504         result += ptr->signed_volume();
00505     return result;
00506 }
00507 
00508 void BSPTreePoly::set_vertex_marks( int value )
00509 {
00510     for( Face* face = faceList; face; face = face->nextPtr )
00511     {
00512         EdgeUse* edge = face->usePtr;
00513         do
00514         {
00515             edge->edgePtr->start()->markVal = value;
00516             edge->edgePtr->end()->markVal   = value;
00517             edge                            = edge->nextPtr;
00518         } while( edge && edge != face->usePtr );
00519     }
00520 }
00521 /*
00522 static void merge_edges( BSPTreePoly::Edge* keep_edge,
00523                          BSPTreePoly::Edge* dead_edge )
00524 {
00525   // edges must share a vertex
00526   BSPTreePoly::Vertex* dead_vtx = keep_edge->common(dead_edge);
00527   assert(dead_vtx);
00528    // vertex may have only two adjacent edges
00529   BSPTreePoly::VertexUse* dead_vtxuse = dead_edge->use(dead_vtx);
00530   assert(dead_vtxuse);
00531   BSPTreePoly::VertexUse* keep_vtxuse = dead_vtxuse->nextPtr;
00532   assert(keep_vtxuse);
00533   assert(keep_vtxuse->edgePtr == keep_edge);
00534   assert(keep_vtxuse->nextPtr == dead_vtxuse);
00535   assert(keep_vtxuse->prevPtr == dead_vtxuse);
00536   assert(dead_vtxuse->prevPtr == keep_vtxuse);
00537 
00538   // kept edge now ends with the kept vertex on the dead edge
00539   keep_vtxuse->set_vertex( dead_edge->other(dead_vtx) );
00540 
00541   // destructors should take care of everything else
00542   // (including removing dead edge from face loops)
00543   delete dead_edge;
00544 }
00545 */
00546 static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx, BSPTreePoly::Edge* into_edge )
00547 {
00548     // split edge, creating new edge
00549     BSPTreePoly::Edge* new_edge = new BSPTreePoly::Edge( new_vtx, into_edge->end() );
00550     into_edge->endPtr->set_vertex( new_vtx );  // This call might delete new_vtx
00551 
00552     // update coedge loops in faces
00553     if( into_edge->forwardPtr )
00554     {
00555         new_edge->forwardPtr = new BSPTreePoly::EdgeUse( new_edge );
00556         new_edge->forwardPtr->insert_after( into_edge->forwardPtr );
00557     }
00558     if( into_edge->reversePtr )
00559     {
00560         new_edge->reversePtr = new BSPTreePoly::EdgeUse( new_edge );
00561         new_edge->reversePtr->insert_before( into_edge->reversePtr );
00562     }
00563 
00564     return new_edge;
00565 }
00566 
00567 static BSPTreePoly::Face* split_face( BSPTreePoly::EdgeUse* start, BSPTreePoly::EdgeUse* end )
00568 {
00569     BSPTreePoly::Face* face = start->facePtr;
00570     assert( face == end->facePtr );
00571     BSPTreePoly::Face* new_face      = new BSPTreePoly::Face;
00572     BSPTreePoly::EdgeUse* keep_start = start->prevPtr;
00573     BSPTreePoly::EdgeUse* keep_end   = end->nextPtr;
00574     for( BSPTreePoly::EdgeUse* ptr = start; ptr != keep_end; ptr = ptr->nextPtr )
00575     {
00576         if( face->usePtr == ptr ) face->usePtr = keep_start;
00577         ptr->facePtr = new_face;
00578     }
00579     new_face->usePtr        = start;
00580     BSPTreePoly::Edge* edge = new BSPTreePoly::Edge( start->start(), end->end() );
00581     edge->forwardPtr        = new BSPTreePoly::EdgeUse( edge );
00582     edge->reversePtr        = new BSPTreePoly::EdgeUse( edge );
00583 
00584     edge->forwardPtr->facePtr = face;
00585     edge->forwardPtr->prevPtr = keep_start;
00586     keep_start->nextPtr       = edge->forwardPtr;
00587     edge->forwardPtr->nextPtr = keep_end;
00588     keep_end->prevPtr         = edge->forwardPtr;
00589 
00590     edge->reversePtr->facePtr = new_face;
00591     edge->reversePtr->nextPtr = start;
00592     start->prevPtr            = edge->reversePtr;
00593     edge->reversePtr->prevPtr = end;
00594     end->nextPtr              = edge->reversePtr;
00595 
00596     return new_face;
00597 }
00598 
00599 bool BSPTreePoly::cut_polyhedron( const CartVect& plane_normal, double plane_coeff )
00600 {
00601     const double EPSILON = 1e-6;  // points this close are considered coincident
00602 
00603     // scale epsilon rather than normalizing normal vector
00604     const double epsilon = EPSILON * ( plane_normal % plane_normal );
00605 
00606     // Classify all points above/below plane and destroy any faces
00607     // that have no vertices below the plane.
00608     const int UNKNOWN = 0;
00609     const int ABOVE   = 1;
00610     const int ON      = 2;
00611     const int BELOW   = 3;
00612     int num_above     = 0;
00613     set_vertex_marks( UNKNOWN );
00614 
00615     // Classify all points above/below plane and
00616     // split any edge that intersect the plane.
00617     for( Face* face = faceList; face; face = face->nextPtr )
00618     {
00619         EdgeUse* edge = face->usePtr;
00620 
00621         do
00622         {
00623             Vertex* start = edge->edgePtr->start();
00624             Vertex* end   = edge->edgePtr->end();
00625 
00626             if( !start->markVal )
00627             {
00628                 double d = plane_normal % *start + plane_coeff;
00629                 if( d * d <= epsilon )
00630                     start->markVal = ON;
00631                 else if( d < 0.0 )
00632                     start->markVal = BELOW;
00633                 else
00634                 {
00635                     start->markVal = ABOVE;
00636                     ++num_above;
00637                 }
00638             }
00639 
00640             if( !end->markVal )
00641             {
00642                 double d = plane_normal % *end + plane_coeff;
00643                 if( d * d <= epsilon )
00644                     end->markVal = ON;
00645                 else if( d < 0.0 )
00646                     end->markVal = BELOW;
00647                 else
00648                 {
00649                     end->markVal = ABOVE;
00650                     ++num_above;
00651                 }
00652             }
00653 
00654             if( ( end->markVal == ABOVE && start->markVal == BELOW ) ||
00655                 ( end->markVal == BELOW && start->markVal == ABOVE ) )
00656             {
00657                 CartVect dir     = *end - *start;
00658                 double t         = -( plane_normal % *start + plane_coeff ) / ( dir % plane_normal );
00659                 Vertex* new_vtx  = new Vertex( *start + t * dir );
00660                 new_vtx->markVal = ON;
00661                 split_edge( new_vtx, edge->edgePtr );  // This call might delete new_vtx
00662                 end = new_vtx;
00663             }
00664 
00665             edge = edge->nextPtr;
00666         } while( edge && edge != face->usePtr );
00667     }
00668 
00669     if( !num_above ) return false;
00670 
00671     // Split faces
00672     for( Face* face = faceList; face; face = face->nextPtr )
00673     {
00674         EdgeUse* edge = face->usePtr;
00675 
00676         EdgeUse *split_start = 0, *split_end = 0, *other_split = 0;
00677         do
00678         {
00679             if( edge->end()->markVal == ON && edge->start()->markVal != ON )
00680             {
00681                 if( !split_start )
00682                     split_start = edge->nextPtr;
00683                 else if( !split_end )
00684                     split_end = edge;
00685                 else
00686                     other_split = edge;
00687             }
00688 
00689             edge = edge->nextPtr;
00690         } while( edge && edge != face->usePtr );
00691 
00692         // If two vertices are on plane (but not every vertex)
00693         // then split the face
00694         if( split_end && !other_split )
00695         {
00696             assert( split_start );
00697             Face* new_face    = split_face( split_start, split_end );
00698             new_face->nextPtr = faceList;
00699             faceList          = new_face;
00700         }
00701     }
00702 
00703     // Destroy all faces that are above the plane
00704     Face** lptr = &faceList;
00705     while( *lptr )
00706     {
00707         EdgeUse* edge   = ( *lptr )->usePtr;
00708         bool some_above = false;
00709         do
00710         {
00711             if( edge->start()->markVal == ABOVE )
00712             {
00713                 some_above = true;
00714                 break;
00715             }
00716             edge = edge->nextPtr;
00717         } while( edge && edge != ( *lptr )->usePtr );
00718 
00719         if( some_above )
00720         {
00721             Face* dead = *lptr;
00722             *lptr      = ( *lptr )->nextPtr;
00723             delete dead;
00724         }
00725         else
00726         {
00727             lptr = &( ( *lptr )->nextPtr );
00728         }
00729     }
00730 
00731     // Construct a new face in the cut plane
00732 
00733     // First find an edge to start at
00734     Edge* edge_ptr = 0;
00735     for( Face* face = faceList; face && !edge_ptr; face = face->nextPtr )
00736     {
00737         EdgeUse* co_edge = face->usePtr;
00738         do
00739         {
00740             if( 0 == co_edge->edgePtr->other( co_edge ) )
00741             {
00742                 edge_ptr = co_edge->edgePtr;
00743                 break;
00744             }
00745             co_edge = co_edge->nextPtr;
00746         } while( co_edge && co_edge != face->usePtr );
00747     }
00748     if( !edge_ptr ) return false;
00749 
00750     // Constuct new face and first CoEdge
00751     faceList = new Face( faceList );
00752     Vertex *next_vtx, *start_vtx;
00753     EdgeUse* prev_coedge;
00754     if( edge_ptr->forwardPtr )
00755     {
00756         next_vtx  = edge_ptr->start();
00757         start_vtx = edge_ptr->end();
00758         assert( !edge_ptr->reversePtr );
00759         prev_coedge = edge_ptr->reversePtr = new EdgeUse( edge_ptr, faceList );
00760     }
00761     else
00762     {
00763         next_vtx    = edge_ptr->end();
00764         start_vtx   = edge_ptr->start();
00765         prev_coedge = edge_ptr->forwardPtr = new EdgeUse( edge_ptr, faceList );
00766     }
00767 
00768     // Construct coedges until loop is closed
00769     while( next_vtx != start_vtx )
00770     {
00771         // find next edge adjacent to vertex with only one adjacent face
00772         VertexUse* this_use = edge_ptr->use( next_vtx );
00773         VertexUse* use      = this_use->nextPtr;
00774         while( use != this_use )
00775         {
00776             if( use->edgePtr->forwardPtr == 0 )
00777             {
00778                 edge_ptr = use->edgePtr;
00779                 assert( edge_ptr->start() == next_vtx );
00780                 next_vtx             = edge_ptr->end();
00781                 edge_ptr->forwardPtr = new EdgeUse( edge_ptr );
00782                 edge_ptr->forwardPtr->insert_after( prev_coedge );
00783                 prev_coedge = edge_ptr->forwardPtr;
00784                 break;
00785             }
00786             else if( use->edgePtr->reversePtr == 0 )
00787             {
00788                 edge_ptr = use->edgePtr;
00789                 assert( edge_ptr->end() == next_vtx );
00790                 next_vtx             = edge_ptr->start();
00791                 edge_ptr->reversePtr = new EdgeUse( edge_ptr );
00792                 edge_ptr->reversePtr->insert_after( prev_coedge );
00793                 prev_coedge = edge_ptr->reversePtr;
00794                 break;
00795             }
00796 
00797             use = use->nextPtr;
00798             assert( use != this_use );  // failed to close loop!
00799         }
00800     }
00801 
00802     return true;
00803 }
00804 
00805 bool BSPTreePoly::is_valid() const
00806 {
00807     std::set< Face* > list_faces;
00808 
00809     int i = 0;
00810     for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
00811     {
00812         if( ++i > 10000 ) return false;
00813         if( !list_faces.insert( ptr ).second ) return false;
00814     }
00815 
00816     std::set< Vertex* > vertices;
00817     for( Face* face = faceList; face; face = face->nextPtr )
00818     {
00819         i               = 0;
00820         EdgeUse* coedge = face->usePtr;
00821         do
00822         {
00823             if( ++i > 10000 ) return false;
00824 
00825             if( coedge->facePtr != face ) return false;
00826 
00827             Edge* edge = coedge->edgePtr;
00828             if( !edge->startPtr || !edge->endPtr ) return false;
00829 
00830             vertices.insert( edge->start() );
00831             vertices.insert( edge->end() );
00832 
00833             EdgeUse* other;
00834             if( edge->forwardPtr == coedge )
00835                 other = edge->reversePtr;
00836             else if( edge->reversePtr != coedge )
00837                 return false;
00838             else
00839                 other = edge->forwardPtr;
00840             if( !other ) return false;
00841             if( list_faces.find( other->facePtr ) == list_faces.end() ) return false;
00842 
00843             EdgeUse* next = coedge->nextPtr;
00844             if( next->prevPtr != coedge ) return false;
00845             if( coedge->end() != next->start() ) return false;
00846 
00847             coedge = next;
00848         } while( coedge != face->usePtr );
00849     }
00850 
00851     for( std::set< Vertex* >::iterator j = vertices.begin(); j != vertices.end(); ++j )
00852     {
00853         Vertex* vtx = *j;
00854 
00855         i              = 0;
00856         VertexUse* use = vtx->usePtr;
00857         do
00858         {
00859             if( ++i > 10000 ) return false;
00860 
00861             if( use->vtxPtr != vtx ) return false;
00862 
00863             Edge* edge = use->edgePtr;
00864             if( !edge ) return false;
00865             if( edge->startPtr != use && edge->endPtr != use ) return false;
00866 
00867             VertexUse* next = use->nextPtr;
00868             if( next->prevPtr != use ) return false;
00869 
00870             use = next;
00871         } while( use != vtx->usePtr );
00872     }
00873 
00874     return true;
00875 }
00876 
00877 bool BSPTreePoly::is_point_contained( const CartVect& point ) const
00878 {
00879     if( !faceList )  // empty (zero-dimension) polyhedron
00880         return false;
00881 
00882     const double EPSILON = 1e-6;
00883     // Test that point is below the plane of each face
00884     // NOTE: This will NOT work for polyhedra w/ concavities
00885     for( Face* face = faceList; face; face = face->nextPtr )
00886     {
00887         Vertex *pt1, *pt2, *pt3;
00888         pt1 = face->usePtr->start();
00889         pt2 = face->usePtr->end();
00890         pt3 = face->usePtr->nextPtr->end();
00891 
00892         if( pt3 == pt1 )  // degenerate
00893             continue;
00894 
00895         CartVect norm = ( *pt3 - *pt2 ) * ( *pt1 - *pt2 );
00896         double coeff  = -( norm % *pt2 );
00897         if( ( norm % point + coeff ) > EPSILON )  // if above plane, with some -epsilon
00898             return false;
00899     }
00900 
00901     return true;
00902 }
00903 
00904 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines