1 #include "moab/CartVect.hpp"
2 #include "moab/BSPTreePoly.hpp"
3 #include <assert.h>
4 #include <stdlib.h>
5 #include <set>
6
7 #undef DEBUG_IDS
8
9 namespace moab {
10
11 struct BSPTreePoly::Vertex : public CartVect {
Vertexmoab::BSPTreePoly::Vertex12 Vertex( const CartVect& v ) : CartVect(v), usePtr(0), markVal(0)
13 #ifdef DEBUG_IDS
14 , id(nextID++)
15 #endif
16 {}
~Vertexmoab::BSPTreePoly::Vertex17 ~Vertex() { assert(!usePtr); }
18 BSPTreePoly::VertexUse* usePtr;
19 int markVal;
20 #ifdef DEBUG_IDS
21 int id;
22 static int nextID;
23 #endif
24 };
25
26 struct BSPTreePoly::VertexUse {
27 VertexUse( Edge* edge, Vertex* vtx );
28 ~VertexUse();
29
30 void set_vertex( BSPTreePoly::Vertex*& vtx_ptr );
31
32 BSPTreePoly::VertexUse *nextPtr, *prevPtr;
33 BSPTreePoly::Vertex* vtxPtr;
34 BSPTreePoly::Edge* edgePtr;
35 };
36
37
38 struct BSPTreePoly::EdgeUse {
39 EdgeUse( Edge* edge );
40 EdgeUse( Edge* edge, Face* face );
41 ~EdgeUse();
42
43 BSPTreePoly::EdgeUse *prevPtr, *nextPtr;
44 BSPTreePoly::Edge* edgePtr;
45 BSPTreePoly::Face* facePtr;
46
47 inline BSPTreePoly::Vertex* start() const;
48 inline BSPTreePoly::Vertex* end() const;
49 int sense() const;
50
51 void insert_after( BSPTreePoly::EdgeUse* prev );
52 void insert_before( BSPTreePoly::EdgeUse* next );
53 };
54
55 struct BSPTreePoly::Edge {
56 BSPTreePoly::VertexUse *startPtr, *endPtr;
57 BSPTreePoly::EdgeUse *forwardPtr, *reversePtr;
58 #ifdef DEBUG_IDS
59 int id;
60 static int nextID;
61 #endif
62
Edgemoab::BSPTreePoly::Edge63 Edge( Vertex* vstart, Vertex* vend )
64 : forwardPtr(0), reversePtr(0)
65 #ifdef DEBUG_IDS
66 , id(nextID++)
67 #endif
68 {
69 startPtr = new VertexUse( this, vstart );
70 endPtr = new VertexUse( this, vend );
71 }
72
73 ~Edge();
74
startmoab::BSPTreePoly::Edge75 BSPTreePoly::Vertex* start() const
76 { return startPtr->vtxPtr; }
endmoab::BSPTreePoly::Edge77 BSPTreePoly::Vertex* end() const
78 { return endPtr->vtxPtr; }
79
forwardmoab::BSPTreePoly::Edge80 BSPTreePoly::Face* forward() const
81 { return forwardPtr ? forwardPtr->facePtr : 0; }
reversemoab::BSPTreePoly::Edge82 BSPTreePoly::Face* reverse() const
83 { return reversePtr ? reversePtr->facePtr : 0; }
84
usemoab::BSPTreePoly::Edge85 BSPTreePoly::VertexUse* use( BSPTreePoly::Vertex* vtx ) const
86 { return (vtx == startPtr->vtxPtr) ? startPtr : (vtx == endPtr->vtxPtr) ? endPtr : 0; }
nextmoab::BSPTreePoly::Edge87 BSPTreePoly::Edge* next( BSPTreePoly::Vertex* about ) const
88 { return use(about)->nextPtr->edgePtr; }
prevmoab::BSPTreePoly::Edge89 BSPTreePoly::Edge* prev( BSPTreePoly::Vertex* about ) const
90 { return use(about)->prevPtr->edgePtr; }
91
usemoab::BSPTreePoly::Edge92 BSPTreePoly::EdgeUse* use( BSPTreePoly::Face* face ) const
93 { return (face == forwardPtr->facePtr) ? forwardPtr : (face == reversePtr->facePtr) ? reversePtr : 0; }
nextmoab::BSPTreePoly::Edge94 BSPTreePoly::Edge* next( BSPTreePoly::Face* about ) const
95 { return use(about)->nextPtr->edgePtr; }
prevmoab::BSPTreePoly::Edge96 BSPTreePoly::Edge* prev( BSPTreePoly::Face* about ) const
97 { return use(about)->prevPtr->edgePtr; }
98
othermoab::BSPTreePoly::Edge99 BSPTreePoly::VertexUse* other( BSPTreePoly::VertexUse* vuse ) const
100 { return vuse == startPtr ? endPtr : vuse == endPtr ? startPtr : 0; }
othermoab::BSPTreePoly::Edge101 BSPTreePoly::EdgeUse* other( BSPTreePoly::EdgeUse* vuse ) const
102 { return vuse == forwardPtr ? reversePtr : vuse == reversePtr ? forwardPtr : 0; }
othermoab::BSPTreePoly::Edge103 BSPTreePoly::Vertex* other( BSPTreePoly::Vertex* vtx ) const
104 { return vtx == startPtr->vtxPtr ? endPtr->vtxPtr :
105 vtx == endPtr->vtxPtr ? startPtr->vtxPtr : 0; }
commonmoab::BSPTreePoly::Edge106 BSPTreePoly::Vertex* common( BSPTreePoly::Edge* eother ) const
107 { return start() == eother->start() || start() == eother->end() ? start() :
108 end() == eother->start() || end() == eother->end() ? end() : 0; }
109
110 int sense( BSPTreePoly::Face* face ) const;
111
112 void remove_from_vertex( BSPTreePoly::Vertex*& vtx_ptr );
113 void remove_from_face( BSPTreePoly::Face*& face_ptr );
114 void add_to_vertex( BSPTreePoly::Vertex* vtx_ptr );
115 };
116
117 struct BSPTreePoly::Face {
Facemoab::BSPTreePoly::Face118 Face(Face* next) : usePtr(0), nextPtr(next)
119 #ifdef DEBUG_IDS
120 , id(nextID++)
121 #endif
122 {}
Facemoab::BSPTreePoly::Face123 Face() : usePtr(0), nextPtr(0)
124 #ifdef DEBUG_IDS
125 , id(nextID++)
126 #endif
127 {}
128 ~Face();
129 BSPTreePoly::EdgeUse* usePtr;
130 BSPTreePoly::Face* nextPtr;
131 #ifdef DEBUG_IDS
132 int id;
133 static int nextID;
134 #endif
135 double signed_volume() const;
136 };
137
138 #ifdef DEBUG_IDS
139 int BSPTreePoly::Vertex::nextID = 1;
140 int BSPTreePoly::Edge::nextID = 1;
141 int BSPTreePoly::Face::nextID = 1;
142 #endif
reset_debug_ids()143 void BSPTreePoly::reset_debug_ids() {
144 #ifdef DEBUG_IDS
145 BSPTreePoly::Vertex::nextID = 1;
146 BSPTreePoly::Edge::nextID = 1;
147 BSPTreePoly::Face::nextID = 1;
148 #endif
149 }
150
151 //static void merge_edges( BSPTreePoly::Edge* keep_edge,
152 // BSPTreePoly::Edge* dead_edge );
153
154 static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx,
155 BSPTreePoly::Edge* into_edge );
156
VertexUse(BSPTreePoly::Edge * edge,BSPTreePoly::Vertex * vtx)157 BSPTreePoly::VertexUse::VertexUse( BSPTreePoly::Edge* edge, BSPTreePoly::Vertex* vtx )
158 : vtxPtr(vtx), edgePtr(edge)
159 {
160 if (!vtx->usePtr) {
161 vtx->usePtr = prevPtr = nextPtr = this;
162 return;
163 }
164
165 nextPtr = vtx->usePtr;
166 prevPtr = nextPtr->prevPtr;
167 assert( prevPtr->nextPtr == nextPtr );
168 nextPtr->prevPtr = this;
169 prevPtr->nextPtr = this;
170 }
171
~VertexUse()172 BSPTreePoly::VertexUse::~VertexUse()
173 {
174 if (nextPtr == this) {
175 assert(prevPtr == this);
176 assert(vtxPtr->usePtr == this);
177 vtxPtr->usePtr = 0;
178 delete vtxPtr;
179 }
180 else if (vtxPtr->usePtr == this)
181 vtxPtr->usePtr = nextPtr;
182
183 nextPtr->prevPtr = prevPtr;
184 prevPtr->nextPtr = nextPtr;
185 nextPtr = prevPtr = 0;
186 }
187
set_vertex(BSPTreePoly::Vertex * & vtx)188 void BSPTreePoly::VertexUse::set_vertex( BSPTreePoly::Vertex*& vtx )
189 {
190 if (vtxPtr) {
191 if (nextPtr == prevPtr) {
192 assert(nextPtr == this);
193 vtxPtr->usePtr = 0;
194 delete vtx;
195 vtx = 0;
196 }
197 else {
198 nextPtr->prevPtr = prevPtr;
199 prevPtr->nextPtr = nextPtr;
200 if (vtxPtr->usePtr == this)
201 vtxPtr->usePtr = nextPtr;
202 }
203 }
204
205 if (vtx) {
206 vtxPtr = vtx;
207 nextPtr = vtxPtr->usePtr->nextPtr;
208 prevPtr = vtxPtr->usePtr;
209 nextPtr->prevPtr = this;
210 vtxPtr->usePtr->nextPtr = this;
211 }
212 }
213
214
EdgeUse(BSPTreePoly::Edge * edge)215 BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge )
216 : prevPtr(0), nextPtr(0), edgePtr(edge), facePtr(0)
217 {}
218
EdgeUse(BSPTreePoly::Edge * edge,BSPTreePoly::Face * face)219 BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge,
220 BSPTreePoly::Face* face )
221 : edgePtr(edge), facePtr(face)
222 {
223 assert(!face->usePtr);
224 face->usePtr = prevPtr = nextPtr = this;
225
226 if (!face->usePtr) {
227 face->usePtr = prevPtr = nextPtr = this;
228 return;
229 }
230
231 nextPtr = face->usePtr;
232 prevPtr = nextPtr->prevPtr;
233 assert( prevPtr->nextPtr == nextPtr );
234 nextPtr->prevPtr = this;
235 prevPtr->nextPtr = this;
236 }
237
insert_after(BSPTreePoly::EdgeUse * prev)238 void BSPTreePoly::EdgeUse::insert_after( BSPTreePoly::EdgeUse* prev )
239 {
240 // shouldn't already be in a face
241 assert(!facePtr);
242 // adjacent edges should share vertices
243 assert( start() == prev->end() );
244
245 facePtr = prev->facePtr;
246 nextPtr = prev->nextPtr;
247 prevPtr = prev;
248 nextPtr->prevPtr = this;
249 prevPtr->nextPtr = this;
250 }
251
insert_before(BSPTreePoly::EdgeUse * next)252 void BSPTreePoly::EdgeUse::insert_before( BSPTreePoly::EdgeUse* next )
253 {
254 // shouldn't already be in a face
255 assert(!facePtr);
256 // adjacent edges should share vertices
257 assert( end() == next->start() );
258
259 facePtr = next->facePtr;
260 prevPtr = next->prevPtr;
261 nextPtr = next;
262 nextPtr->prevPtr = this;
263 prevPtr->nextPtr = this;
264 }
265
266
~EdgeUse()267 BSPTreePoly::EdgeUse::~EdgeUse()
268 {
269 if (facePtr->usePtr == this)
270 facePtr->usePtr = (nextPtr == this) ? 0 : nextPtr;
271
272 if (edgePtr->forwardPtr == this)
273 edgePtr->forwardPtr = 0;
274 if (edgePtr->reversePtr == this)
275 edgePtr->reversePtr = 0;
276
277 if (!edgePtr->forwardPtr && !edgePtr->reversePtr)
278 delete edgePtr;
279
280 nextPtr->prevPtr = prevPtr;
281 prevPtr->nextPtr = nextPtr;
282 nextPtr = prevPtr = 0;
283 }
284
sense() const285 int BSPTreePoly::EdgeUse::sense() const
286 {
287 if (edgePtr->forwardPtr == this)
288 return 1;
289 else if (edgePtr->reversePtr == this)
290 return -1;
291 else
292 return 0;
293 }
294
start() const295 BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::start() const
296 {
297 if (edgePtr->forwardPtr == this)
298 return edgePtr->start();
299 else if (edgePtr->reversePtr == this)
300 return edgePtr->end();
301 else
302 return 0;
303 }
304
end() const305 BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::end() const
306 {
307 if (edgePtr->forwardPtr == this)
308 return edgePtr->end();
309 else if (edgePtr->reversePtr == this)
310 return edgePtr->start();
311 else
312 return 0;
313 }
314
~Edge()315 BSPTreePoly::Edge::~Edge()
316 {
317 delete startPtr;
318 delete endPtr;
319 delete forwardPtr;
320 delete reversePtr;
321 }
322
sense(BSPTreePoly::Face * face) const323 int BSPTreePoly::Edge::sense( BSPTreePoly::Face* face ) const
324 {
325 if (forwardPtr && forwardPtr->facePtr == face)
326 return 1;
327 else if (reversePtr && reversePtr->facePtr == face)
328 return -1;
329 else
330 return 0;
331 }
332
~Face()333 BSPTreePoly::Face::~Face()
334 {
335 BSPTreePoly::EdgeUse* nextEdgeUsePtr = usePtr;
336 while (nextEdgeUsePtr) {
337 delete nextEdgeUsePtr; // This is tricky: ~EdgeUse() might change the value of usePtr
338 if (usePtr && usePtr != nextEdgeUsePtr)
339 nextEdgeUsePtr = usePtr;
340 else
341 nextEdgeUsePtr = 0;
342 }
343 usePtr = 0;
344 }
345
clear()346 void BSPTreePoly::clear() {
347 while (faceList) {
348 Face* face = faceList;
349 faceList = faceList->nextPtr;
350 delete face;
351 }
352 }
353
set(const CartVect hex_corners[8])354 ErrorCode BSPTreePoly::set( const CartVect hex_corners[8] )
355 {
356 clear();
357
358 Vertex* vertices[8];
359 for (int i = 0; i < 8; ++i)
360 vertices[i] = new Vertex(hex_corners[i]);
361
362 Edge* edges[12];
363 #ifdef DEBUG_IDS
364 int start_id = Edge::nextID;
365 #endif
366 for (int i = 0; i < 4; ++i) {
367 int j= (i+1)%4;
368 edges[i ] = new Edge( vertices[i ], vertices[j ] );
369 edges[i+4] = new Edge( vertices[i ], vertices[i+4] );
370 edges[i+8] = new Edge( vertices[i+4], vertices[j+4] );
371 }
372 #ifdef DEBUG_IDS
373 for (int i = 0; i < 12; ++i)
374 edges[i]->id = start_id++;
375 #endif
376
377 static const int face_conn[6][4] = { { 0, 5, -8, -4 },
378 { 1, 6, -9, -5 },
379 { 2, 7, -10, -6 },
380 { 3, 4, -11, -7 },
381 { -3, -2, -1,-12 },
382 { 8, 9, 10, 11 } };
383 for (int i = 0; i < 6; ++i) {
384 faceList = new Face(faceList);
385 EdgeUse* prev = 0;
386 for (int j = 0; j < 4; ++j) {
387 int e = face_conn[i][j];
388 if (e < 0) {
389 e = (-e) % 12;
390 assert(!edges[e]->reversePtr);
391 if (!prev) {
392 edges[e]->reversePtr = new EdgeUse( edges[e], faceList );
393 }
394 else {
395 edges[e]->reversePtr = new EdgeUse( edges[e] );
396 edges[e]->reversePtr->insert_after( prev );
397 }
398 prev = edges[e]->reversePtr;
399 }
400 else {
401 assert(!edges[e]->forwardPtr);
402 if (!prev) {
403 edges[e]->forwardPtr = new EdgeUse( edges[e], faceList );
404 }
405 else {
406 edges[e]->forwardPtr = new EdgeUse( edges[e] );
407 edges[e]->forwardPtr->insert_after( prev );
408 }
409 prev = edges[e]->forwardPtr;
410 }
411 }
412 }
413
414 return MB_SUCCESS;
415 }
416
get_faces(std::vector<const Face * > & face_list) const417 void BSPTreePoly::get_faces( std::vector<const Face*>& face_list ) const
418 {
419 face_list.clear();
420 for (Face* face = faceList; face; face = face->nextPtr)
421 face_list.push_back( face );
422 }
423
get_vertices(const Face * face,std::vector<CartVect> & vertices) const424 void BSPTreePoly::get_vertices( const Face* face,
425 std::vector<CartVect>& vertices ) const
426 {
427 vertices.clear();
428 if (!face || !face->usePtr)
429 return;
430
431 EdgeUse* coedge = face->usePtr;
432 do {
433 vertices.push_back( *coedge->end() );
434 coedge = coedge->nextPtr;
435 } while (coedge != face->usePtr);
436 }
437
signed_volume() const438 double BSPTreePoly::Face::signed_volume() const
439 {
440 CartVect sum(0.0);
441 const CartVect* base = usePtr->start();
442 CartVect d1 = (*usePtr->end() - *base);
443 for (EdgeUse* coedge = usePtr->nextPtr; coedge != usePtr; coedge = coedge->nextPtr) {
444 CartVect d2 = (*coedge->end() - *base);
445 sum += d1 * d2;
446 d1 = d2;
447 }
448 return (1.0/6.0) * (sum % *base);
449 }
450
volume() const451 double BSPTreePoly::volume() const
452 {
453 double result = 0;
454 for (Face* ptr = faceList; ptr; ptr = ptr->nextPtr)
455 result += ptr->signed_volume();
456 return result;
457 }
458
set_vertex_marks(int value)459 void BSPTreePoly::set_vertex_marks( int value )
460 {
461 for (Face* face = faceList; face; face = face->nextPtr) {
462 EdgeUse* edge = face->usePtr;
463 do {
464 edge->edgePtr->start()->markVal = value;
465 edge->edgePtr->end()->markVal = value;
466 edge = edge->nextPtr;
467 } while (edge && edge != face->usePtr);
468 }
469 }
470 /*
471 static void merge_edges( BSPTreePoly::Edge* keep_edge,
472 BSPTreePoly::Edge* dead_edge )
473 {
474 // edges must share a vertex
475 BSPTreePoly::Vertex* dead_vtx = keep_edge->common(dead_edge);
476 assert(dead_vtx);
477 // vertex may have only two adjacent edges
478 BSPTreePoly::VertexUse* dead_vtxuse = dead_edge->use(dead_vtx);
479 assert(dead_vtxuse);
480 BSPTreePoly::VertexUse* keep_vtxuse = dead_vtxuse->nextPtr;
481 assert(keep_vtxuse);
482 assert(keep_vtxuse->edgePtr == keep_edge);
483 assert(keep_vtxuse->nextPtr == dead_vtxuse);
484 assert(keep_vtxuse->prevPtr == dead_vtxuse);
485 assert(dead_vtxuse->prevPtr == keep_vtxuse);
486
487 // kept edge now ends with the kept vertex on the dead edge
488 keep_vtxuse->set_vertex( dead_edge->other(dead_vtx) );
489
490 // destructors should take care of everything else
491 // (including removing dead edge from face loops)
492 delete dead_edge;
493 }
494 */
split_edge(BSPTreePoly::Vertex * & new_vtx,BSPTreePoly::Edge * into_edge)495 static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx,
496 BSPTreePoly::Edge* into_edge )
497 {
498 // split edge, creating new edge
499 BSPTreePoly::Edge* new_edge = new BSPTreePoly::Edge( new_vtx, into_edge->end() );
500 into_edge->endPtr->set_vertex( new_vtx ); // This call might delete new_vtx
501
502 // update coedge loops in faces
503 if (into_edge->forwardPtr) {
504 new_edge->forwardPtr = new BSPTreePoly::EdgeUse( new_edge );
505 new_edge->forwardPtr->insert_after( into_edge->forwardPtr );
506 }
507 if (into_edge->reversePtr) {
508 new_edge->reversePtr = new BSPTreePoly::EdgeUse( new_edge );
509 new_edge->reversePtr->insert_before( into_edge->reversePtr );
510 }
511
512 return new_edge;
513 }
514
split_face(BSPTreePoly::EdgeUse * start,BSPTreePoly::EdgeUse * end)515 static BSPTreePoly::Face* split_face( BSPTreePoly::EdgeUse* start,
516 BSPTreePoly::EdgeUse* end )
517 {
518 BSPTreePoly::Face* face = start->facePtr;
519 assert(face == end->facePtr);
520 BSPTreePoly::Face* new_face = new BSPTreePoly::Face;
521 BSPTreePoly::EdgeUse* keep_start = start->prevPtr;
522 BSPTreePoly::EdgeUse* keep_end = end->nextPtr;
523 for (BSPTreePoly::EdgeUse* ptr = start; ptr != keep_end; ptr = ptr->nextPtr) {
524 if (face->usePtr == ptr)
525 face->usePtr = keep_start;
526 ptr->facePtr = new_face;
527 }
528 new_face->usePtr = start;
529 BSPTreePoly::Edge* edge = new BSPTreePoly::Edge( start->start(), end->end() );
530 edge->forwardPtr = new BSPTreePoly::EdgeUse( edge );
531 edge->reversePtr = new BSPTreePoly::EdgeUse( edge );
532
533 edge->forwardPtr->facePtr = face;
534 edge->forwardPtr->prevPtr = keep_start;
535 keep_start->nextPtr = edge->forwardPtr;
536 edge->forwardPtr->nextPtr = keep_end;
537 keep_end->prevPtr = edge->forwardPtr;
538
539 edge->reversePtr->facePtr = new_face;
540 edge->reversePtr->nextPtr = start;
541 start->prevPtr = edge->reversePtr;
542 edge->reversePtr->prevPtr = end;
543 end->nextPtr = edge->reversePtr;
544
545 return new_face;
546 }
547
548
cut_polyhedron(const CartVect & plane_normal,double plane_coeff)549 bool BSPTreePoly::cut_polyhedron( const CartVect& plane_normal,
550 double plane_coeff )
551 {
552 const double EPSILON = 1e-6; // points this close are considered coincident
553
554 // scale epsilon rather than normalizing normal vector
555 const double epsilon = EPSILON * (plane_normal % plane_normal);
556
557 // Classify all points above/below plane and destroy any faces
558 // that have no vertices below the plane.
559 const int UNKNOWN = 0;
560 const int ABOVE = 1;
561 const int ON = 2;
562 const int BELOW = 3;
563 int num_above = 0;
564 set_vertex_marks( UNKNOWN );
565
566 // Classify all points above/below plane and
567 // split any edge that intersect the plane.
568 for (Face* face = faceList; face; face = face->nextPtr) {
569 EdgeUse* edge = face->usePtr;
570
571 do {
572 Vertex* start = edge->edgePtr->start();
573 Vertex* end = edge->edgePtr->end();
574
575 if (!start->markVal) {
576 double d = plane_normal % *start + plane_coeff;
577 if (d*d <= epsilon)
578 start->markVal = ON;
579 else if (d < 0.0)
580 start->markVal = BELOW;
581 else {
582 start->markVal = ABOVE;
583 ++num_above;
584 }
585 }
586
587 if (!end->markVal) {
588 double d = plane_normal % *end + plane_coeff;
589 if (d*d <= epsilon)
590 end->markVal = ON;
591 else if (d < 0.0)
592 end->markVal = BELOW;
593 else {
594 end->markVal = ABOVE;
595 ++num_above;
596 }
597 }
598
599 if ((end->markVal == ABOVE && start->markVal == BELOW) ||
600 (end->markVal == BELOW && start->markVal == ABOVE)) {
601 CartVect dir = *end - *start;
602 double t = -(plane_normal % *start + plane_coeff) / (dir % plane_normal);
603 Vertex* new_vtx = new Vertex( *start + t*dir );
604 new_vtx->markVal = ON;
605 split_edge( new_vtx, edge->edgePtr ); // This call might delete new_vtx
606 end = new_vtx;
607 }
608
609 edge = edge->nextPtr;
610 } while (edge && edge != face->usePtr);
611 }
612
613 if (!num_above)
614 return false;
615
616 // Split faces
617 for (Face* face = faceList; face; face = face->nextPtr) {
618 EdgeUse* edge = face->usePtr;
619
620 EdgeUse *split_start = 0, *split_end = 0, *other_split = 0;
621 do {
622 if (edge->end()->markVal == ON && edge->start()->markVal != ON) {
623 if (!split_start)
624 split_start = edge->nextPtr;
625 else if (!split_end)
626 split_end = edge;
627 else
628 other_split = edge;
629 }
630
631 edge = edge->nextPtr;
632 } while (edge && edge != face->usePtr);
633
634 // If two vertices are on plane (but not every vertex)
635 // then split the face
636 if (split_end && !other_split) {
637 assert(split_start);
638 Face* new_face = split_face( split_start, split_end );
639 new_face->nextPtr = faceList;
640 faceList = new_face;
641 }
642 }
643
644 // Destroy all faces that are above the plane
645 Face** lptr = &faceList;
646 while (*lptr) {
647 EdgeUse* edge = (*lptr)->usePtr;
648 bool some_above = false;
649 do {
650 if (edge->start()->markVal == ABOVE) {
651 some_above = true;
652 break;
653 }
654 edge = edge->nextPtr;
655 } while (edge && edge != (*lptr)->usePtr);
656
657 if (some_above) {
658 Face* dead = *lptr;
659 *lptr = (*lptr)->nextPtr;
660 delete dead;
661 }
662 else {
663 lptr = &((*lptr)->nextPtr);
664 }
665 }
666
667 // Construct a new face in the cut plane
668
669 // First find an edge to start at
670 Edge* edge_ptr = 0;
671 for (Face* face = faceList; face && !edge_ptr; face = face->nextPtr) {
672 EdgeUse* co_edge = face->usePtr;
673 do {
674 if (0 == co_edge->edgePtr->other(co_edge)) {
675 edge_ptr = co_edge->edgePtr;
676 break;
677 }
678 co_edge = co_edge->nextPtr;
679 } while (co_edge && co_edge != face->usePtr);
680 }
681 if (!edge_ptr) return false;
682
683 // Constuct new face and first CoEdge
684 faceList = new Face( faceList );
685 Vertex *next_vtx, *start_vtx;
686 EdgeUse* prev_coedge;
687 if (edge_ptr->forwardPtr) {
688 next_vtx = edge_ptr->start();
689 start_vtx = edge_ptr->end();
690 assert(!edge_ptr->reversePtr);
691 prev_coedge = edge_ptr->reversePtr = new EdgeUse( edge_ptr, faceList );
692 }
693 else {
694 next_vtx = edge_ptr->end();
695 start_vtx = edge_ptr->start();
696 prev_coedge = edge_ptr->forwardPtr = new EdgeUse( edge_ptr, faceList );
697 }
698
699 // Construct coedges until loop is closed
700 while (next_vtx != start_vtx) {
701 // find next edge adjacent to vertex with only one adjacent face
702 VertexUse* this_use = edge_ptr->use( next_vtx );
703 VertexUse* use = this_use->nextPtr;
704 while (use != this_use) {
705 if (use->edgePtr->forwardPtr == 0) {
706 edge_ptr = use->edgePtr;
707 assert(edge_ptr->start() == next_vtx);
708 next_vtx = edge_ptr->end();
709 edge_ptr->forwardPtr = new EdgeUse( edge_ptr );
710 edge_ptr->forwardPtr->insert_after( prev_coedge );
711 prev_coedge = edge_ptr->forwardPtr;
712 break;
713 }
714 else if (use->edgePtr->reversePtr == 0) {
715 edge_ptr = use->edgePtr;
716 assert(edge_ptr->end() == next_vtx);
717 next_vtx = edge_ptr->start();
718 edge_ptr->reversePtr = new EdgeUse( edge_ptr );
719 edge_ptr->reversePtr->insert_after( prev_coedge );
720 prev_coedge = edge_ptr->reversePtr;
721 break;
722 }
723
724 use = use->nextPtr;
725 assert( use != this_use ); // failed to close loop!
726 }
727 }
728
729 return true;
730 }
731
is_valid() const732 bool BSPTreePoly::is_valid() const
733 {
734 std::set<Face*> list_faces;
735
736 int i = 0;
737 for (Face* ptr = faceList; ptr; ptr = ptr->nextPtr) {
738 if (++i > 10000)
739 return false;
740 if (!list_faces.insert(ptr).second)
741 return false;
742 }
743
744 std::set<Vertex*> vertices;
745 for (Face* face = faceList; face; face = face->nextPtr) {
746 i = 0;
747 EdgeUse* coedge = face->usePtr;
748 do {
749 if (++i > 10000)
750 return false;
751
752 if (coedge->facePtr != face)
753 return false;
754
755 Edge* edge = coedge->edgePtr;
756 if (!edge->startPtr || !edge->endPtr)
757 return false;
758
759 vertices.insert( edge->start() );
760 vertices.insert( edge->end() );
761
762 EdgeUse* other;
763 if (edge->forwardPtr == coedge)
764 other = edge->reversePtr;
765 else if (edge->reversePtr != coedge)
766 return false;
767 else
768 other = edge->forwardPtr;
769 if (!other)
770 return false;
771 if (list_faces.find( other->facePtr ) == list_faces.end())
772 return false;
773
774 EdgeUse* next = coedge->nextPtr;
775 if (next->prevPtr != coedge)
776 return false;
777 if (coedge->end() != next->start())
778 return false;
779
780 coedge = next;
781 } while (coedge != face->usePtr);
782 }
783
784 for (std::set<Vertex*>::iterator j = vertices.begin(); j != vertices.end(); ++j) {
785 Vertex* vtx = *j;
786
787 i = 0;
788 VertexUse* use = vtx->usePtr;
789 do {
790 if (++i > 10000)
791 return false;
792
793 if (use->vtxPtr != vtx)
794 return false;
795
796 Edge* edge = use->edgePtr;
797 if (!edge)
798 return false;
799 if (edge->startPtr != use && edge->endPtr != use)
800 return false;
801
802 VertexUse* next = use->nextPtr;
803 if (next->prevPtr != use)
804 return false;
805
806 use = next;
807 } while (use != vtx->usePtr);
808 }
809
810 return true;
811 }
812
is_point_contained(const CartVect & point) const813 bool BSPTreePoly::is_point_contained( const CartVect& point ) const
814 {
815 if (!faceList) // empty (zero-dimension) polyhedron
816 return false;
817
818 const double EPSILON = 1e-6;
819 // Test that point is below the plane of each face
820 // NOTE: This will NOT work for polyhedra w/ concavities
821 for (Face* face = faceList; face; face = face->nextPtr) {
822 Vertex *pt1, *pt2, *pt3;
823 pt1 = face->usePtr->start();
824 pt2 = face->usePtr->end();
825 pt3 = face->usePtr->nextPtr->end();
826
827 if (pt3 == pt1) // degenerate
828 continue;
829
830 CartVect norm = (*pt3 - *pt2) * (*pt1 - *pt2);
831 double coeff = -(norm % *pt2);
832 if ( (norm % point + coeff) > EPSILON) // if above plane, with some -epsilon
833 return false;
834 }
835
836 return true;
837 }
838
839 } // namespace moab
840