MOAB: Mesh Oriented datABase  (version 5.2.1)
BSPTreePoly.cpp
Go to the documentation of this file.
00001 #include "moab/CartVect.hpp"
00002 #include "moab/BSPTreePoly.hpp"
00003 #include <assert.h>
00004 #include <stdlib.h>
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()
00145                    ? start()
00146                    : end() == eother->start() || end() == eother->end() ? end() : 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 ) { edges[e]->reversePtr = new EdgeUse( edges[e], faceList ); }
00436                 else
00437                 {
00438                     edges[e]->reversePtr = new EdgeUse( edges[e] );
00439                     edges[e]->reversePtr->insert_after( prev );
00440                 }
00441                 prev = edges[e]->reversePtr;
00442             }
00443             else
00444             {
00445                 assert( !edges[e]->forwardPtr );
00446                 if( !prev ) { edges[e]->forwardPtr = new EdgeUse( edges[e], faceList ); }
00447                 else
00448                 {
00449                     edges[e]->forwardPtr = new EdgeUse( edges[e] );
00450                     edges[e]->forwardPtr->insert_after( prev );
00451                 }
00452                 prev = edges[e]->forwardPtr;
00453             }
00454         }
00455     }
00456 
00457     return MB_SUCCESS;
00458 }
00459 
00460 void BSPTreePoly::get_faces( std::vector< const Face* >& face_list ) const
00461 {
00462     face_list.clear();
00463     for( Face* face = faceList; face; face = face->nextPtr )
00464         face_list.push_back( face );
00465 }
00466 
00467 void BSPTreePoly::get_vertices( const Face* face, std::vector< CartVect >& vertices ) const
00468 {
00469     vertices.clear();
00470     if( !face || !face->usePtr ) return;
00471 
00472     EdgeUse* coedge = face->usePtr;
00473     do
00474     {
00475         vertices.push_back( *coedge->end() );
00476         coedge = coedge->nextPtr;
00477     } while( coedge != face->usePtr );
00478 }
00479 
00480 double BSPTreePoly::Face::signed_volume() const
00481 {
00482     CartVect sum( 0.0 );
00483     const CartVect* base = usePtr->start();
00484     CartVect d1          = ( *usePtr->end() - *base );
00485     for( EdgeUse* coedge = usePtr->nextPtr; coedge != usePtr; coedge = coedge->nextPtr )
00486     {
00487         CartVect d2 = ( *coedge->end() - *base );
00488         sum += d1 * d2;
00489         d1 = d2;
00490     }
00491     return ( 1.0 / 6.0 ) * ( sum % *base );
00492 }
00493 
00494 double BSPTreePoly::volume() const
00495 {
00496     double result = 0;
00497     for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
00498         result += ptr->signed_volume();
00499     return result;
00500 }
00501 
00502 void BSPTreePoly::set_vertex_marks( int value )
00503 {
00504     for( Face* face = faceList; face; face = face->nextPtr )
00505     {
00506         EdgeUse* edge = face->usePtr;
00507         do
00508         {
00509             edge->edgePtr->start()->markVal = value;
00510             edge->edgePtr->end()->markVal   = value;
00511             edge                            = edge->nextPtr;
00512         } while( edge && edge != face->usePtr );
00513     }
00514 }
00515 /*
00516 static void merge_edges( BSPTreePoly::Edge* keep_edge,
00517                          BSPTreePoly::Edge* dead_edge )
00518 {
00519   // edges must share a vertex
00520   BSPTreePoly::Vertex* dead_vtx = keep_edge->common(dead_edge);
00521   assert(dead_vtx);
00522    // vertex may have only two adjacent edges
00523   BSPTreePoly::VertexUse* dead_vtxuse = dead_edge->use(dead_vtx);
00524   assert(dead_vtxuse);
00525   BSPTreePoly::VertexUse* keep_vtxuse = dead_vtxuse->nextPtr;
00526   assert(keep_vtxuse);
00527   assert(keep_vtxuse->edgePtr == keep_edge);
00528   assert(keep_vtxuse->nextPtr == dead_vtxuse);
00529   assert(keep_vtxuse->prevPtr == dead_vtxuse);
00530   assert(dead_vtxuse->prevPtr == keep_vtxuse);
00531 
00532   // kept edge now ends with the kept vertex on the dead edge
00533   keep_vtxuse->set_vertex( dead_edge->other(dead_vtx) );
00534 
00535   // destructors should take care of everything else
00536   // (including removing dead edge from face loops)
00537   delete dead_edge;
00538 }
00539 */
00540 static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx, BSPTreePoly::Edge* into_edge )
00541 {
00542     // split edge, creating new edge
00543     BSPTreePoly::Edge* new_edge = new BSPTreePoly::Edge( new_vtx, into_edge->end() );
00544     into_edge->endPtr->set_vertex( new_vtx );  // This call might delete new_vtx
00545 
00546     // update coedge loops in faces
00547     if( into_edge->forwardPtr )
00548     {
00549         new_edge->forwardPtr = new BSPTreePoly::EdgeUse( new_edge );
00550         new_edge->forwardPtr->insert_after( into_edge->forwardPtr );
00551     }
00552     if( into_edge->reversePtr )
00553     {
00554         new_edge->reversePtr = new BSPTreePoly::EdgeUse( new_edge );
00555         new_edge->reversePtr->insert_before( into_edge->reversePtr );
00556     }
00557 
00558     return new_edge;
00559 }
00560 
00561 static BSPTreePoly::Face* split_face( BSPTreePoly::EdgeUse* start, BSPTreePoly::EdgeUse* end )
00562 {
00563     BSPTreePoly::Face* face = start->facePtr;
00564     assert( face == end->facePtr );
00565     BSPTreePoly::Face* new_face      = new BSPTreePoly::Face;
00566     BSPTreePoly::EdgeUse* keep_start = start->prevPtr;
00567     BSPTreePoly::EdgeUse* keep_end   = end->nextPtr;
00568     for( BSPTreePoly::EdgeUse* ptr = start; ptr != keep_end; ptr = ptr->nextPtr )
00569     {
00570         if( face->usePtr == ptr ) face->usePtr = keep_start;
00571         ptr->facePtr = new_face;
00572     }
00573     new_face->usePtr        = start;
00574     BSPTreePoly::Edge* edge = new BSPTreePoly::Edge( start->start(), end->end() );
00575     edge->forwardPtr        = new BSPTreePoly::EdgeUse( edge );
00576     edge->reversePtr        = new BSPTreePoly::EdgeUse( edge );
00577 
00578     edge->forwardPtr->facePtr = face;
00579     edge->forwardPtr->prevPtr = keep_start;
00580     keep_start->nextPtr       = edge->forwardPtr;
00581     edge->forwardPtr->nextPtr = keep_end;
00582     keep_end->prevPtr         = edge->forwardPtr;
00583 
00584     edge->reversePtr->facePtr = new_face;
00585     edge->reversePtr->nextPtr = start;
00586     start->prevPtr            = edge->reversePtr;
00587     edge->reversePtr->prevPtr = end;
00588     end->nextPtr              = edge->reversePtr;
00589 
00590     return new_face;
00591 }
00592 
00593 bool BSPTreePoly::cut_polyhedron( const CartVect& plane_normal, double plane_coeff )
00594 {
00595     const double EPSILON = 1e-6;  // points this close are considered coincident
00596 
00597     // scale epsilon rather than normalizing normal vector
00598     const double epsilon = EPSILON * ( plane_normal % plane_normal );
00599 
00600     // Classify all points above/below plane and destroy any faces
00601     // that have no vertices below the plane.
00602     const int UNKNOWN = 0;
00603     const int ABOVE   = 1;
00604     const int ON      = 2;
00605     const int BELOW   = 3;
00606     int num_above     = 0;
00607     set_vertex_marks( UNKNOWN );
00608 
00609     // Classify all points above/below plane and
00610     // split any edge that intersect the plane.
00611     for( Face* face = faceList; face; face = face->nextPtr )
00612     {
00613         EdgeUse* edge = face->usePtr;
00614 
00615         do
00616         {
00617             Vertex* start = edge->edgePtr->start();
00618             Vertex* end   = edge->edgePtr->end();
00619 
00620             if( !start->markVal )
00621             {
00622                 double d = plane_normal % *start + plane_coeff;
00623                 if( d * d <= epsilon )
00624                     start->markVal = ON;
00625                 else if( d < 0.0 )
00626                     start->markVal = BELOW;
00627                 else
00628                 {
00629                     start->markVal = ABOVE;
00630                     ++num_above;
00631                 }
00632             }
00633 
00634             if( !end->markVal )
00635             {
00636                 double d = plane_normal % *end + plane_coeff;
00637                 if( d * d <= epsilon )
00638                     end->markVal = ON;
00639                 else if( d < 0.0 )
00640                     end->markVal = BELOW;
00641                 else
00642                 {
00643                     end->markVal = ABOVE;
00644                     ++num_above;
00645                 }
00646             }
00647 
00648             if( ( end->markVal == ABOVE && start->markVal == BELOW ) ||
00649                 ( end->markVal == BELOW && start->markVal == ABOVE ) )
00650             {
00651                 CartVect dir     = *end - *start;
00652                 double t         = -( plane_normal % *start + plane_coeff ) / ( dir % plane_normal );
00653                 Vertex* new_vtx  = new Vertex( *start + t * dir );
00654                 new_vtx->markVal = ON;
00655                 split_edge( new_vtx, edge->edgePtr );  // This call might delete new_vtx
00656                 end = new_vtx;
00657             }
00658 
00659             edge = edge->nextPtr;
00660         } while( edge && edge != face->usePtr );
00661     }
00662 
00663     if( !num_above ) return false;
00664 
00665     // Split faces
00666     for( Face* face = faceList; face; face = face->nextPtr )
00667     {
00668         EdgeUse* edge = face->usePtr;
00669 
00670         EdgeUse *split_start = 0, *split_end = 0, *other_split = 0;
00671         do
00672         {
00673             if( edge->end()->markVal == ON && edge->start()->markVal != ON )
00674             {
00675                 if( !split_start )
00676                     split_start = edge->nextPtr;
00677                 else if( !split_end )
00678                     split_end = edge;
00679                 else
00680                     other_split = edge;
00681             }
00682 
00683             edge = edge->nextPtr;
00684         } while( edge && edge != face->usePtr );
00685 
00686         // If two vertices are on plane (but not every vertex)
00687         // then split the face
00688         if( split_end && !other_split )
00689         {
00690             assert( split_start );
00691             Face* new_face    = split_face( split_start, split_end );
00692             new_face->nextPtr = faceList;
00693             faceList          = new_face;
00694         }
00695     }
00696 
00697     // Destroy all faces that are above the plane
00698     Face** lptr = &faceList;
00699     while( *lptr )
00700     {
00701         EdgeUse* edge   = ( *lptr )->usePtr;
00702         bool some_above = false;
00703         do
00704         {
00705             if( edge->start()->markVal == ABOVE )
00706             {
00707                 some_above = true;
00708                 break;
00709             }
00710             edge = edge->nextPtr;
00711         } while( edge && edge != ( *lptr )->usePtr );
00712 
00713         if( some_above )
00714         {
00715             Face* dead = *lptr;
00716             *lptr      = ( *lptr )->nextPtr;
00717             delete dead;
00718         }
00719         else
00720         {
00721             lptr = &( ( *lptr )->nextPtr );
00722         }
00723     }
00724 
00725     // Construct a new face in the cut plane
00726 
00727     // First find an edge to start at
00728     Edge* edge_ptr = 0;
00729     for( Face* face = faceList; face && !edge_ptr; face = face->nextPtr )
00730     {
00731         EdgeUse* co_edge = face->usePtr;
00732         do
00733         {
00734             if( 0 == co_edge->edgePtr->other( co_edge ) )
00735             {
00736                 edge_ptr = co_edge->edgePtr;
00737                 break;
00738             }
00739             co_edge = co_edge->nextPtr;
00740         } while( co_edge && co_edge != face->usePtr );
00741     }
00742     if( !edge_ptr ) return false;
00743 
00744     // Constuct new face and first CoEdge
00745     faceList = new Face( faceList );
00746     Vertex *next_vtx, *start_vtx;
00747     EdgeUse* prev_coedge;
00748     if( edge_ptr->forwardPtr )
00749     {
00750         next_vtx  = edge_ptr->start();
00751         start_vtx = edge_ptr->end();
00752         assert( !edge_ptr->reversePtr );
00753         prev_coedge = edge_ptr->reversePtr = new EdgeUse( edge_ptr, faceList );
00754     }
00755     else
00756     {
00757         next_vtx    = edge_ptr->end();
00758         start_vtx   = edge_ptr->start();
00759         prev_coedge = edge_ptr->forwardPtr = new EdgeUse( edge_ptr, faceList );
00760     }
00761 
00762     // Construct coedges until loop is closed
00763     while( next_vtx != start_vtx )
00764     {
00765         // find next edge adjacent to vertex with only one adjacent face
00766         VertexUse* this_use = edge_ptr->use( next_vtx );
00767         VertexUse* use      = this_use->nextPtr;
00768         while( use != this_use )
00769         {
00770             if( use->edgePtr->forwardPtr == 0 )
00771             {
00772                 edge_ptr = use->edgePtr;
00773                 assert( edge_ptr->start() == next_vtx );
00774                 next_vtx             = edge_ptr->end();
00775                 edge_ptr->forwardPtr = new EdgeUse( edge_ptr );
00776                 edge_ptr->forwardPtr->insert_after( prev_coedge );
00777                 prev_coedge = edge_ptr->forwardPtr;
00778                 break;
00779             }
00780             else if( use->edgePtr->reversePtr == 0 )
00781             {
00782                 edge_ptr = use->edgePtr;
00783                 assert( edge_ptr->end() == next_vtx );
00784                 next_vtx             = edge_ptr->start();
00785                 edge_ptr->reversePtr = new EdgeUse( edge_ptr );
00786                 edge_ptr->reversePtr->insert_after( prev_coedge );
00787                 prev_coedge = edge_ptr->reversePtr;
00788                 break;
00789             }
00790 
00791             use = use->nextPtr;
00792             assert( use != this_use );  // failed to close loop!
00793         }
00794     }
00795 
00796     return true;
00797 }
00798 
00799 bool BSPTreePoly::is_valid() const
00800 {
00801     std::set< Face* > list_faces;
00802 
00803     int i = 0;
00804     for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
00805     {
00806         if( ++i > 10000 ) return false;
00807         if( !list_faces.insert( ptr ).second ) return false;
00808     }
00809 
00810     std::set< Vertex* > vertices;
00811     for( Face* face = faceList; face; face = face->nextPtr )
00812     {
00813         i               = 0;
00814         EdgeUse* coedge = face->usePtr;
00815         do
00816         {
00817             if( ++i > 10000 ) return false;
00818 
00819             if( coedge->facePtr != face ) return false;
00820 
00821             Edge* edge = coedge->edgePtr;
00822             if( !edge->startPtr || !edge->endPtr ) return false;
00823 
00824             vertices.insert( edge->start() );
00825             vertices.insert( edge->end() );
00826 
00827             EdgeUse* other;
00828             if( edge->forwardPtr == coedge )
00829                 other = edge->reversePtr;
00830             else if( edge->reversePtr != coedge )
00831                 return false;
00832             else
00833                 other = edge->forwardPtr;
00834             if( !other ) return false;
00835             if( list_faces.find( other->facePtr ) == list_faces.end() ) return false;
00836 
00837             EdgeUse* next = coedge->nextPtr;
00838             if( next->prevPtr != coedge ) return false;
00839             if( coedge->end() != next->start() ) return false;
00840 
00841             coedge = next;
00842         } while( coedge != face->usePtr );
00843     }
00844 
00845     for( std::set< Vertex* >::iterator j = vertices.begin(); j != vertices.end(); ++j )
00846     {
00847         Vertex* vtx = *j;
00848 
00849         i              = 0;
00850         VertexUse* use = vtx->usePtr;
00851         do
00852         {
00853             if( ++i > 10000 ) return false;
00854 
00855             if( use->vtxPtr != vtx ) return false;
00856 
00857             Edge* edge = use->edgePtr;
00858             if( !edge ) return false;
00859             if( edge->startPtr != use && edge->endPtr != use ) return false;
00860 
00861             VertexUse* next = use->nextPtr;
00862             if( next->prevPtr != use ) return false;
00863 
00864             use = next;
00865         } while( use != vtx->usePtr );
00866     }
00867 
00868     return true;
00869 }
00870 
00871 bool BSPTreePoly::is_point_contained( const CartVect& point ) const
00872 {
00873     if( !faceList )  // empty (zero-dimension) polyhedron
00874         return false;
00875 
00876     const double EPSILON = 1e-6;
00877     // Test that point is below the plane of each face
00878     // NOTE: This will NOT work for polyhedra w/ concavities
00879     for( Face* face = faceList; face; face = face->nextPtr )
00880     {
00881         Vertex *pt1, *pt2, *pt3;
00882         pt1 = face->usePtr->start();
00883         pt2 = face->usePtr->end();
00884         pt3 = face->usePtr->nextPtr->end();
00885 
00886         if( pt3 == pt1 )  // degenerate
00887             continue;
00888 
00889         CartVect norm = ( *pt3 - *pt2 ) * ( *pt1 - *pt2 );
00890         double coeff  = -( norm % *pt2 );
00891         if( ( norm % point + coeff ) > EPSILON )  // if above plane, with some -epsilon
00892             return false;
00893     }
00894 
00895     return true;
00896 }
00897 
00898 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines