MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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