[MOAB-dev] r3168 - MOAB/trunk
kraftche at cae.wisc.edu
kraftche at cae.wisc.edu
Fri Sep 25 10:34:25 CDT 2009
Author: kraftche
Date: 2009-09-25 10:34:25 -0500 (Fri, 25 Sep 2009)
New Revision: 3168
Added:
MOAB/trunk/BSPTreePoly.cpp
MOAB/trunk/BSPTreePoly.hpp
Modified:
MOAB/trunk/MBAdaptiveKDTree.cpp
MOAB/trunk/MBAdaptiveKDTree.hpp
MOAB/trunk/MBBSPTree.cpp
MOAB/trunk/MBBSPTree.hpp
MOAB/trunk/MBGeomUtil.cpp
MOAB/trunk/MBGeomUtil.hpp
MOAB/trunk/Makefile.am
MOAB/trunk/adaptive_kd_tree_tests.cpp
MOAB/trunk/bsp_tree_test.cpp
Log:
o Add function to test for intersection between a ray and a kD-tree leaf,
and return the points at which it enters and exits the leaf.
o Add function to test for intersection between a ray and a BSP-tree leaf,
and return the points at which it enters and exits the leaf.
o Add fast version of BSP-tree ray/leaf intersect for tree where all
leaves are boxes.
o Add a class for calculating the polyhedron representing the boundary
of a leaf of a general BSP-tree.
o Add function to calculate the volume of a leaf of a general BSP-tree
by first calculating the bounding polyhedron.
o Add unit tests for the above
Added: MOAB/trunk/BSPTreePoly.cpp
===================================================================
--- MOAB/trunk/BSPTreePoly.cpp (rev 0)
+++ MOAB/trunk/BSPTreePoly.cpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -0,0 +1,823 @@
+#include "MBCartVect.hpp"
+#include "BSPTreePoly.hpp"
+#include <assert.h>
+#include <stdlib.h>
+#include <set>
+
+#undef DEBUG_IDS
+
+struct BSPTreePoly::Vertex : public MBCartVect {
+ Vertex( const MBCartVect& v ) : MBCartVect(v), usePtr(0)
+#ifdef DEBUG_IDS
+ , id(nextID++)
+#endif
+ {}
+ ~Vertex() { assert(!usePtr); }
+ BSPTreePoly::VertexUse* usePtr;
+ int markVal;
+#ifdef DEBUG_IDS
+ int id;
+ static int nextID;
+#endif
+};
+
+struct BSPTreePoly::VertexUse {
+ VertexUse( Edge* edge, Vertex* vtx );
+ ~VertexUse();
+
+ void set_vertex( BSPTreePoly::Vertex* vtx_ptr );
+
+ BSPTreePoly::VertexUse *nextPtr, *prevPtr;
+ BSPTreePoly::Vertex* vtxPtr;
+ BSPTreePoly::Edge* edgePtr;
+};
+
+
+struct BSPTreePoly::EdgeUse {
+ EdgeUse( Edge* edge );
+ EdgeUse( Edge* edge, Face* face );
+ ~EdgeUse();
+
+ BSPTreePoly::EdgeUse *prevPtr, *nextPtr;
+ BSPTreePoly::Edge* edgePtr;
+ BSPTreePoly::Face* facePtr;
+
+ inline BSPTreePoly::Vertex* start() const;
+ inline BSPTreePoly::Vertex* end() const;
+ int sense() const;
+
+ void insert_after( BSPTreePoly::EdgeUse* prev );
+ void insert_before( BSPTreePoly::EdgeUse* next );
+};
+
+struct BSPTreePoly::Edge {
+ BSPTreePoly::VertexUse *startPtr, *endPtr;
+ BSPTreePoly::EdgeUse *forwardPtr, *reversePtr;
+#ifdef DEBUG_IDS
+ int id;
+ static int nextID;
+#endif
+
+ Edge( Vertex* start, Vertex* end )
+ : forwardPtr(0), reversePtr(0)
+#ifdef DEBUG_IDS
+ , id(nextID++)
+#endif
+ {
+ startPtr = new VertexUse( this, start );
+ endPtr = new VertexUse( this, end );
+ }
+
+ ~Edge();
+
+ BSPTreePoly::Vertex* start() const
+ { return startPtr->vtxPtr; }
+ BSPTreePoly::Vertex* end() const
+ { return endPtr->vtxPtr; }
+
+ BSPTreePoly::Face* forward() const
+ { return forwardPtr ? forwardPtr->facePtr : 0; }
+ BSPTreePoly::Face* reverse() const
+ { return reversePtr ? reversePtr->facePtr : 0; }
+
+ BSPTreePoly::VertexUse* use( BSPTreePoly::Vertex* vtx ) const
+ { return (vtx == startPtr->vtxPtr) ? startPtr : (vtx == endPtr->vtxPtr) ? endPtr : 0; }
+ BSPTreePoly::Edge* next( BSPTreePoly::Vertex* about ) const
+ { return use(about)->nextPtr->edgePtr; }
+ BSPTreePoly::Edge* prev( BSPTreePoly::Vertex* about ) const
+ { return use(about)->prevPtr->edgePtr; }
+
+ BSPTreePoly::EdgeUse* use( BSPTreePoly::Face* face ) const
+ { return (face == forwardPtr->facePtr) ? forwardPtr : (face == reversePtr->facePtr) ? reversePtr : 0; }
+ BSPTreePoly::Edge* next( BSPTreePoly::Face* about ) const
+ { return use(about)->nextPtr->edgePtr; }
+ BSPTreePoly::Edge* prev( BSPTreePoly::Face* about ) const
+ { return use(about)->prevPtr->edgePtr; }
+
+ BSPTreePoly::VertexUse* other( BSPTreePoly::VertexUse* use ) const
+ { return use == startPtr ? endPtr : use == endPtr ? startPtr : 0; }
+ BSPTreePoly::EdgeUse* other( BSPTreePoly::EdgeUse* use ) const
+ { return use == forwardPtr ? reversePtr : use == reversePtr ? forwardPtr : 0; }
+ BSPTreePoly::Vertex* other( BSPTreePoly::Vertex* vtx ) const
+ { return vtx == startPtr->vtxPtr ? endPtr->vtxPtr :
+ vtx == endPtr->vtxPtr ? startPtr->vtxPtr : 0; }
+ BSPTreePoly::Vertex* common( BSPTreePoly::Edge* other ) const
+ { return start() == other->start() || start() == other->end() ? start() :
+ end() == other->start() || end() == other->end() ? end() : 0; }
+
+ int sense( BSPTreePoly::Face* face ) const;
+
+ void remove_from_vertex( BSPTreePoly::Vertex*& vtx_ptr );
+ void remove_from_face( BSPTreePoly::Face*& face_ptr );
+ void add_to_vertex( BSPTreePoly::Vertex* vtx_ptr );
+};
+
+struct BSPTreePoly::Face {
+ Face(Face* next) : usePtr(0), nextPtr(next)
+#ifdef DEBUG_IDS
+ , id(nextID++)
+#endif
+ {}
+ Face() : usePtr(0), nextPtr(0)
+#ifdef DEBUG_IDS
+ , id(nextID++)
+#endif
+ {}
+ ~Face();
+ BSPTreePoly::EdgeUse* usePtr;
+ BSPTreePoly::Face* nextPtr;
+#ifdef DEBUG_IDS
+ int id;
+ static int nextID;
+#endif
+ double signed_volume() const;
+};
+
+#ifdef DEBUG_IDS
+int BSPTreePoly::Vertex::nextID = 1;
+int BSPTreePoly::Edge::nextID = 1;
+int BSPTreePoly::Face::nextID = 1;
+#endif
+void BSPTreePoly::reset_debug_ids() {
+#ifdef DEBUG_IDS
+ BSPTreePoly::Vertex::nextID = 1;
+ BSPTreePoly::Edge::nextID = 1;
+ BSPTreePoly::Face::nextID = 1;
+#endif
+}
+
+static void merge_edges( BSPTreePoly::Edge* keep_edge,
+ BSPTreePoly::Edge* dead_edge );
+
+static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex* new_vtx,
+ BSPTreePoly::Edge* into_edge );
+
+BSPTreePoly::VertexUse::VertexUse( BSPTreePoly::Edge* edge, BSPTreePoly::Vertex* vtx )
+ : vtxPtr(vtx), edgePtr(edge)
+{
+ if (!vtx->usePtr) {
+ vtx->usePtr = prevPtr = nextPtr = this;
+ return;
+ }
+
+ nextPtr = vtx->usePtr;
+ prevPtr = nextPtr->prevPtr;
+ assert( prevPtr->nextPtr == nextPtr );
+ nextPtr->prevPtr = this;
+ prevPtr->nextPtr = this;
+}
+
+BSPTreePoly::VertexUse::~VertexUse()
+{
+ if (nextPtr == this) {
+ assert(prevPtr == this);
+ assert(vtxPtr->usePtr == this);
+ vtxPtr->usePtr = 0;
+ delete vtxPtr;
+ }
+ else if (vtxPtr->usePtr == this)
+ vtxPtr->usePtr = nextPtr;
+
+ nextPtr->prevPtr = prevPtr;
+ prevPtr->nextPtr = nextPtr;
+ nextPtr = prevPtr = 0;
+}
+
+void BSPTreePoly::VertexUse::set_vertex( BSPTreePoly::Vertex* vtx )
+{
+ if (vtxPtr) {
+ if (nextPtr == prevPtr) {
+ assert(nextPtr == this);
+ vtxPtr->usePtr = 0;
+ delete vtx;
+ }
+ else {
+ nextPtr->prevPtr = prevPtr;
+ prevPtr->nextPtr = nextPtr;
+ if (vtxPtr->usePtr == this)
+ vtxPtr->usePtr = nextPtr;
+ }
+ }
+
+ vtxPtr = vtx;
+ nextPtr = vtxPtr->usePtr->nextPtr;
+ prevPtr = vtxPtr->usePtr;
+ nextPtr->prevPtr = this;
+ vtxPtr->usePtr->nextPtr = this;
+}
+
+
+BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge )
+ : prevPtr(0), nextPtr(0), edgePtr(edge), facePtr(0)
+{}
+
+BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge,
+ BSPTreePoly::Face* face )
+ : edgePtr(edge), facePtr(face)
+{
+ assert(!face->usePtr);
+ face->usePtr = prevPtr = nextPtr = this;
+
+ if (!face->usePtr) {
+ face->usePtr = prevPtr = nextPtr = this;
+ return;
+ }
+
+ nextPtr = face->usePtr;
+ prevPtr = nextPtr->prevPtr;
+ assert( prevPtr->nextPtr == nextPtr );
+ nextPtr->prevPtr = this;
+ prevPtr->nextPtr = this;
+}
+
+void BSPTreePoly::EdgeUse::insert_after( BSPTreePoly::EdgeUse* prev )
+{
+ // shouldn't aready be in a face
+ assert(!facePtr);
+ // adjacent edges should share vertices
+ assert( start() == prev->end() );
+
+ facePtr = prev->facePtr;
+ nextPtr = prev->nextPtr;
+ prevPtr = prev;
+ nextPtr->prevPtr = this;
+ prevPtr->nextPtr = this;
+}
+
+void BSPTreePoly::EdgeUse::insert_before( BSPTreePoly::EdgeUse* next )
+{
+ // shouldn't aready be in a face
+ assert(!facePtr);
+ // adjacent edges should share vertices
+ assert( end() == next->start() );
+
+ facePtr = next->facePtr;
+ prevPtr = next->prevPtr;
+ nextPtr = next;
+ nextPtr->prevPtr = this;
+ prevPtr->nextPtr = this;
+}
+
+
+BSPTreePoly::EdgeUse::~EdgeUse()
+{
+ if (facePtr->usePtr == this)
+ facePtr->usePtr = (nextPtr == this) ? 0 : nextPtr;
+
+ if (edgePtr->forwardPtr == this)
+ edgePtr->forwardPtr = 0;
+ if (edgePtr->reversePtr == this)
+ edgePtr->reversePtr = 0;
+
+ if (!edgePtr->forwardPtr && !edgePtr->reversePtr)
+ delete edgePtr;
+
+ nextPtr->prevPtr = prevPtr;
+ prevPtr->nextPtr = nextPtr;
+ nextPtr = prevPtr = 0;
+}
+
+int BSPTreePoly::EdgeUse::sense() const
+{
+ if (edgePtr->forwardPtr == this)
+ return 1;
+ else if (edgePtr->reversePtr == this)
+ return -1;
+ else
+ return 0;
+};
+
+BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::start() const
+{
+ if (edgePtr->forwardPtr == this)
+ return edgePtr->start();
+ else if (edgePtr->reversePtr == this)
+ return edgePtr->end();
+ else
+ return 0;
+}
+
+BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::end() const
+{
+ if (edgePtr->forwardPtr == this)
+ return edgePtr->end();
+ else if (edgePtr->reversePtr == this)
+ return edgePtr->start();
+ else
+ return 0;
+}
+
+BSPTreePoly::Edge::~Edge()
+{
+ delete startPtr;
+ delete endPtr;
+ delete forwardPtr;
+ delete reversePtr;
+}
+
+int BSPTreePoly::Edge::sense( BSPTreePoly::Face* face ) const
+{
+ if (forwardPtr && forwardPtr->facePtr == face)
+ return 1;
+ else if (reversePtr && reversePtr->facePtr == face)
+ return -1;
+ else
+ return 0;
+}
+
+BSPTreePoly::Face::~Face()
+{
+ while (usePtr)
+ delete usePtr;
+}
+
+void BSPTreePoly::clear() {
+ while (faceList) {
+ Face* face = faceList;
+ faceList = faceList->nextPtr;
+ delete face;
+ }
+}
+
+MBErrorCode BSPTreePoly::set( const MBCartVect hex_corners[8] )
+{
+ clear();
+
+ Vertex* vertices[8];
+ for (int i = 0; i < 8; ++i)
+ vertices[i] = new Vertex(hex_corners[i]);
+
+ Edge* edges[12];
+#ifdef DEBUG_IDS
+ int start_id = Edge::nextID;
+#endif
+ for (int i = 0; i < 4; ++i) {
+ int j= (i+1)%4;
+ edges[i ] = new Edge( vertices[i ], vertices[j ] );
+ edges[i+4] = new Edge( vertices[i ], vertices[i+4] );
+ edges[i+8] = new Edge( vertices[i+4], vertices[j+4] );
+ }
+#ifdef DEBUG_IDS
+ for (int i = 0; i < 12; ++i)
+ edges[i]->id = start_id++;
+#endif
+
+ static const int face_conn[6][4] = { { 0, 5, -8, -4 },
+ { 1, 6, -9, -5 },
+ { 2, 7, -10, -6 },
+ { 3, 4, -11, -7 },
+ { -3, -2, -1,-12 },
+ { 8, 9, 10, 11 } };
+ for (int i = 0; i < 6; ++i) {
+ faceList = new Face(faceList);
+ EdgeUse* prev = 0;
+ for (int j = 0; j < 4; ++j) {
+ int e = face_conn[i][j];
+ if (e < 0) {
+ e = (-e) % 12;
+ assert(!edges[e]->reversePtr);
+ if (!prev) {
+ edges[e]->reversePtr = new EdgeUse( edges[e], faceList );
+ }
+ else {
+ edges[e]->reversePtr = new EdgeUse( edges[e] );
+ edges[e]->reversePtr->insert_after( prev );
+ }
+ prev = edges[e]->reversePtr;
+ }
+ else {
+ assert(!edges[e]->forwardPtr);
+ if (!prev) {
+ edges[e]->forwardPtr = new EdgeUse( edges[e], faceList );
+ }
+ else {
+ edges[e]->forwardPtr = new EdgeUse( edges[e] );
+ edges[e]->forwardPtr->insert_after( prev );
+ }
+ prev = edges[e]->forwardPtr;
+ }
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+void BSPTreePoly::get_faces( std::vector<const Face*>& face_list ) const
+{
+ face_list.clear();
+ for (Face* face = faceList; face; face = face->nextPtr)
+ face_list.push_back( face );
+}
+
+void BSPTreePoly::get_vertices( const Face* face,
+ std::vector<MBCartVect>& vertices ) const
+{
+ vertices.clear();
+ if (!face || !face->usePtr)
+ return;
+
+ EdgeUse* coedge = face->usePtr;
+ do {
+ vertices.push_back( *coedge->end() );
+ coedge = coedge->nextPtr;
+ } while (coedge != face->usePtr);
+}
+
+double BSPTreePoly::Face::signed_volume() const
+{
+ MBCartVect sum(0.0);
+ const MBCartVect* base = usePtr->start();
+ MBCartVect d1 = (*usePtr->end() - *base);
+ for (EdgeUse* coedge = usePtr->nextPtr; coedge != usePtr; coedge = coedge->nextPtr) {
+ MBCartVect d2 = (*coedge->end() - *base);
+ sum += d1 * d2;
+ d1 = d2;
+ }
+ return (1.0/6.0) * (sum % *base);
+}
+
+double BSPTreePoly::volume() const
+{
+ double result = 0;
+ for (Face* ptr = faceList; ptr; ptr = ptr->nextPtr)
+ result += ptr->signed_volume();
+ return result;
+}
+
+void BSPTreePoly::set_vertex_marks( int value )
+{
+ for (Face* face = faceList; face; face = face->nextPtr) {
+ EdgeUse* edge = face->usePtr;
+ do {
+ edge->edgePtr->start()->markVal = value;
+ edge->edgePtr->end()->markVal = value;
+ edge = edge->nextPtr;
+ } while (edge && edge != face->usePtr);
+ }
+}
+
+static void merge_edges( BSPTreePoly::Edge* keep_edge,
+ BSPTreePoly::Edge* dead_edge )
+{
+ // edges must share a vertex
+ BSPTreePoly::Vertex* dead_vtx = keep_edge->common(dead_edge);
+ assert(dead_vtx);
+ // vertex may have only two adjacent edges
+ BSPTreePoly::VertexUse* dead_vtxuse = dead_edge->use(dead_vtx);
+ assert(dead_vtxuse);
+ BSPTreePoly::VertexUse* keep_vtxuse = dead_vtxuse->nextPtr;
+ assert(keep_vtxuse);
+ assert(keep_vtxuse->edgePtr == keep_edge);
+ assert(keep_vtxuse->nextPtr == dead_vtxuse);
+ assert(keep_vtxuse->prevPtr == dead_vtxuse);
+ assert(dead_vtxuse->prevPtr == keep_vtxuse);
+
+ // kept edge now ends with the kept vertex on the dead edge
+ keep_vtxuse->set_vertex( dead_edge->other(dead_vtx) );
+
+ // destructors should take care of everything else
+ // (including removing dead edge from face loops)
+ delete dead_edge;
+}
+
+static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex* new_vtx,
+ BSPTreePoly::Edge* into_edge )
+{
+ // split edge, creating new edge
+ BSPTreePoly::Edge* new_edge = new BSPTreePoly::Edge( new_vtx, into_edge->end() );
+ into_edge->endPtr->set_vertex( new_vtx );
+
+ // update coedge loops in faces
+ if (into_edge->forwardPtr) {
+ new_edge->forwardPtr = new BSPTreePoly::EdgeUse( new_edge );
+ new_edge->forwardPtr->insert_after( into_edge->forwardPtr );
+ }
+ if (into_edge->reversePtr) {
+ new_edge->reversePtr = new BSPTreePoly::EdgeUse( new_edge );
+ new_edge->reversePtr->insert_before( into_edge->reversePtr );
+ }
+
+ return new_edge;
+}
+
+static BSPTreePoly::Face* split_face( BSPTreePoly::EdgeUse* start,
+ BSPTreePoly::EdgeUse* end )
+{
+ BSPTreePoly::Face* face = start->facePtr;
+ assert(face == end->facePtr);
+ BSPTreePoly::Face* new_face = new BSPTreePoly::Face;
+ BSPTreePoly::EdgeUse* keep_start = start->prevPtr;
+ BSPTreePoly::EdgeUse* keep_end = end->nextPtr;
+ for (BSPTreePoly::EdgeUse* ptr = start; ptr != keep_end; ptr = ptr->nextPtr) {
+ if (face->usePtr == ptr)
+ face->usePtr = keep_start;
+ ptr->facePtr = new_face;
+ }
+ new_face->usePtr = start;
+ BSPTreePoly::Edge* edge = new BSPTreePoly::Edge( start->start(), end->end() );
+ edge->forwardPtr = new BSPTreePoly::EdgeUse( edge );
+ edge->reversePtr = new BSPTreePoly::EdgeUse( edge );
+
+ edge->forwardPtr->facePtr = face;
+ edge->forwardPtr->prevPtr = keep_start;
+ keep_start->nextPtr = edge->forwardPtr;
+ edge->forwardPtr->nextPtr = keep_end;
+ keep_end->prevPtr = edge->forwardPtr;
+
+ edge->reversePtr->facePtr = new_face;
+ edge->reversePtr->nextPtr = start;
+ start->prevPtr = edge->reversePtr;
+ edge->reversePtr->prevPtr = end;
+ end->nextPtr = edge->reversePtr;
+
+ return new_face;
+}
+
+
+bool BSPTreePoly::cut_polyhedron( const MBCartVect& plane_normal,
+ double plane_coeff )
+{
+ const double EPSILON = 1e-6; // points this close are considered coincident
+
+ // scale epsilon rather than normalizing normal vector
+ const double epsilon = EPSILON * (plane_normal % plane_normal);
+
+ // Classify all points above/below plane and destroy any faces
+ // that have no vertices below the plane.
+ const int UNKNOWN = 0;
+ const int ABOVE = 1;
+ const int ON = 2;
+ const int BELOW = 3;
+ int num_above = 0;
+ set_vertex_marks( UNKNOWN );
+
+ // Classify all points above/below plane and
+ // split any edge that intersect the plane.
+ for (Face* face = faceList; face; face = face->nextPtr) {
+ EdgeUse* edge = face->usePtr;
+
+ do {
+ Vertex* start = edge->edgePtr->start();
+ Vertex* end = edge->edgePtr->end();
+
+ if (!start->markVal) {
+ double d = plane_normal % *start + plane_coeff;
+ if (d*d <= epsilon)
+ start->markVal = ON;
+ else if (d < 0.0)
+ start->markVal = BELOW;
+ else {
+ start->markVal = ABOVE;
+ ++num_above;
+ }
+ }
+
+ if (!end->markVal) {
+ double d = plane_normal % *end + plane_coeff;
+ if (d*d <= epsilon)
+ end->markVal = ON;
+ else if (d < 0.0)
+ end->markVal = BELOW;
+ else {
+ end->markVal = ABOVE;
+ ++num_above;
+ }
+ }
+
+ if ((end->markVal == ABOVE && start->markVal == BELOW) ||
+ (end->markVal == BELOW && start->markVal == ABOVE)) {
+ MBCartVect dir = *end - *start;
+ double t = -(plane_normal % *start + plane_coeff) / (dir % plane_normal);
+ Vertex* new_vtx = new Vertex( *start + t*dir );
+ new_vtx->markVal = ON;
+ split_edge( new_vtx, edge->edgePtr );
+ end = new_vtx;
+ }
+
+ edge = edge->nextPtr;
+ } while (edge && edge != face->usePtr);
+ }
+
+ if (!num_above)
+ return false;
+
+ // Split faces
+ for (Face* face = faceList; face; face = face->nextPtr) {
+ EdgeUse* edge = face->usePtr;
+
+ EdgeUse *split_start = 0, *split_end = 0, *other_split = 0;
+ do {
+ if (edge->end()->markVal == ON && edge->start()->markVal != ON) {
+ if (!split_start)
+ split_start = edge->nextPtr;
+ else if (!split_end)
+ split_end = edge;
+ else
+ other_split = edge;
+ }
+
+ edge = edge->nextPtr;
+ } while (edge && edge != face->usePtr);
+
+ // If two vertices are on plane (but not every vertex)
+ // then split the face
+ if (split_end && !other_split) {
+ assert(split_start);
+ Face* new_face = split_face( split_start, split_end );
+ new_face->nextPtr = faceList;
+ faceList = new_face;
+ }
+ }
+
+ // Destroy all faces that are above the plane
+ Face** lptr = &faceList;
+ while (*lptr) {
+ EdgeUse* edge = (*lptr)->usePtr;
+ bool some_above = false;
+ do {
+ if (edge->start()->markVal == ABOVE) {
+ some_above = true;
+ break;
+ }
+ edge = edge->nextPtr;
+ } while (edge && edge != (*lptr)->usePtr);
+
+ if (some_above) {
+ Face* dead = *lptr;
+ *lptr = (*lptr)->nextPtr;
+ delete dead;
+ }
+ else {
+ lptr = &((*lptr)->nextPtr);
+ }
+ }
+
+ // Construct a new face in the cut plane
+
+ // First find an edge to start at
+ Edge* edge_ptr = 0;
+ for (Face* face = faceList; face && !edge_ptr; face = face->nextPtr) {
+ EdgeUse* co_edge = face->usePtr;
+ do {
+ if (0 == co_edge->edgePtr->other(co_edge)) {
+ edge_ptr = co_edge->edgePtr;
+ break;
+ }
+ co_edge = co_edge->nextPtr;
+ } while (co_edge && co_edge != face->usePtr);
+ }
+
+ // Constuct new face and first CoEdge
+ faceList = new Face( faceList );
+ Vertex *next_vtx, *start_vtx;
+ EdgeUse* prev_coedge;
+ if (edge_ptr->forwardPtr) {
+ next_vtx = edge_ptr->start();
+ start_vtx = edge_ptr->end();
+ assert(!edge_ptr->reversePtr);
+ prev_coedge = edge_ptr->reversePtr = new EdgeUse( edge_ptr, faceList );
+ }
+ else {
+ next_vtx = edge_ptr->end();
+ start_vtx = edge_ptr->start();
+ prev_coedge = edge_ptr->forwardPtr = new EdgeUse( edge_ptr, faceList );
+ }
+
+ // Construct coedges until loop is closed
+ while (next_vtx != start_vtx) {
+ // find next edge adjacent to vertex with only one adjacent face
+ VertexUse* this_use = edge_ptr->use( next_vtx );
+ VertexUse* use = this_use->nextPtr;
+ while (use != this_use) {
+ if (use->edgePtr->forwardPtr == 0) {
+ edge_ptr = use->edgePtr;
+ assert(edge_ptr->start() == next_vtx);
+ next_vtx = edge_ptr->end();
+ edge_ptr->forwardPtr = new EdgeUse( edge_ptr );
+ edge_ptr->forwardPtr->insert_after( prev_coedge );
+ prev_coedge = edge_ptr->forwardPtr;
+ break;
+ }
+ else if (use->edgePtr->reversePtr == 0) {
+ edge_ptr = use->edgePtr;
+ assert(edge_ptr->end() == next_vtx);
+ next_vtx = edge_ptr->start();
+ edge_ptr->reversePtr = new EdgeUse( edge_ptr );
+ edge_ptr->reversePtr->insert_after( prev_coedge );
+ prev_coedge = edge_ptr->reversePtr;
+ break;
+ }
+
+ use = use->nextPtr;
+ assert( use != this_use ); // failed to close loop!
+ }
+ }
+
+ return true;
+}
+
+bool BSPTreePoly::is_valid() const
+{
+ std::set<Face*> list_faces;
+
+ int i = 0;
+ for (Face* ptr = faceList; ptr; ptr = ptr->nextPtr) {
+ if (++i > 10000)
+ return false;
+ if (!list_faces.insert(ptr).second)
+ return false;
+ }
+
+ std::set<Vertex*> vertices;
+ for (Face* face = faceList; face; face = face->nextPtr) {
+ i = 0;
+ EdgeUse* coedge = face->usePtr;
+ do {
+ if (++i > 10000)
+ return false;
+
+ if (coedge->facePtr != face)
+ return false;
+
+ Edge* edge = coedge->edgePtr;
+ if (!edge->startPtr || !edge->endPtr)
+ return false;
+
+ vertices.insert( edge->start() );
+ vertices.insert( edge->end() );
+
+ EdgeUse* other;
+ if (edge->forwardPtr == coedge)
+ other = edge->reversePtr;
+ else if (edge->reversePtr != coedge)
+ return false;
+ else
+ other = edge->forwardPtr;
+ if (!other)
+ return false;
+ if (list_faces.find( other->facePtr ) == list_faces.end())
+ return false;
+
+ EdgeUse* next = coedge->nextPtr;
+ if (next->prevPtr != coedge)
+ return false;
+ if (coedge->end() != next->start())
+ return false;
+
+ coedge = next;
+ } while (coedge != face->usePtr);
+ }
+
+ for (std::set<Vertex*>::iterator j = vertices.begin(); j != vertices.end(); ++j) {
+ Vertex* vtx = *j;
+
+ i = 0;
+ VertexUse* use = vtx->usePtr;
+ do {
+ if (++i > 10000)
+ return false;
+
+ if (use->vtxPtr != vtx)
+ return false;
+
+ Edge* edge = use->edgePtr;
+ if (!edge)
+ return false;
+ if (edge->startPtr != use && edge->endPtr != use)
+ return false;
+
+ VertexUse* next = use->nextPtr;
+ if (next->prevPtr != use)
+ return false;
+
+ use = next;
+ } while (use != vtx->usePtr);
+ }
+
+ return true;
+}
+
+bool BSPTreePoly::is_point_contained( const MBCartVect& point ) const
+{
+ if (!faceList) // empty (zero-dimension) polyhedron
+ return false;
+
+ // Test that point is below the plane of each face
+ // NOTE: This will NOT work for polyhedra w/ concavities
+ for (Face* face = faceList; face; face = face->nextPtr) {
+ Vertex *pt1, *pt2, *pt3;
+ pt1 = face->usePtr->start();
+ pt2 = face->usePtr->end();
+ pt3 = face->usePtr->nextPtr->end();
+
+ if (pt3 == pt1) // degenerate
+ continue;
+
+ MBCartVect norm = (*pt3 - *pt2) * (*pt1 - *pt2);
+ double coeff = -(norm % *pt2);
+ if (norm % point > -coeff) // if above plane
+ return false;
+ }
+
+ return true;
+}
Added: MOAB/trunk/BSPTreePoly.hpp
===================================================================
--- MOAB/trunk/BSPTreePoly.hpp (rev 0)
+++ MOAB/trunk/BSPTreePoly.hpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -0,0 +1,75 @@
+#ifndef MB_BSP_TREE_POLY_HPP
+#define MB_BSP_TREE_POLY_HPP
+
+#include "MBTypes.h"
+#include <vector>
+
+class MBCartVect;
+
+/**\brief Convex polyhedron
+ *
+ * This class is used to represent the convex polyhedron that bounds
+ * a node in a general plane-based BSP-tree.
+ */
+class BSPTreePoly
+{
+ public:
+ struct Vertex;
+ struct VertexUse;
+ struct Edge;
+ struct EdgeUse;
+ struct Face;
+ private:
+ Face* faceList;
+
+ void set_vertex_marks( int value );
+
+ BSPTreePoly( const BSPTreePoly& copy ); // not implemented
+ BSPTreePoly& operator=( const BSPTreePoly& copy ); // not implemented
+
+ public:
+
+ /**\brief Initialize as a planar-faced hexahedron
+ *\param hex_corners Corner coordinates for a hexahedron, in Exodus/Patran order
+ */
+ BSPTreePoly( const MBCartVect hex_corners[8] ) : faceList(0) { set(hex_corners); }
+ BSPTreePoly( ) : faceList(0) { }
+ ~BSPTreePoly() { clear(); }
+
+ /**\brief Initialize as a planar-faced hexahedron
+ *\param hex_corners Corner coordinates for a hexahedron, in Exodus/Patran order
+ */
+ MBErrorCode set( const MBCartVect hex_corners[8] );
+ void clear();
+
+ /**\brief Get handles for faces */
+ void get_faces( std::vector<const Face*>& face_list ) const;
+ /**\brief Get corner coordinates for a face */
+ void get_vertices( const Face* face,
+ std::vector<MBCartVect>& vertices ) const;
+
+ /** Intersect a plane with a polyhedron, retaining
+ * the portion of the polyhedron below the plane.
+ * This will fail if polyhedron is not convex.
+ */
+ bool cut_polyhedron( const MBCartVect& plane_normal,
+ double plane_coeff );
+
+ /** Test if a point is contained in the polyhedron.
+ *
+ *\NOTE algorithm assumes *convex* polyhedron.
+ */
+ bool is_point_contained( const MBCartVect& point ) const;
+
+ //! Assumes planar faces
+ double volume() const;
+
+ // Check that data structure is consistent and represents
+ // a closed polyhedron
+ bool is_valid() const;
+
+ /** For debugging, does nothing unless debug feature is enabled */
+ static void reset_debug_ids();
+};
+
+#endif
Modified: MOAB/trunk/MBAdaptiveKDTree.cpp
===================================================================
--- MOAB/trunk/MBAdaptiveKDTree.cpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/MBAdaptiveKDTree.cpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -752,6 +752,18 @@
return childVect[0] == handle();
}
+bool MBAdaptiveKDTreeIter::intersect_ray( const double ray_point[3],
+ const double ray_vect[3],
+ double& t_enter,
+ double& t_exit ) const
+{
+ return MBGeomUtil::ray_box_intersect( MBCartVect(box_min()),
+ MBCartVect(box_max()),
+ MBCartVect(ray_point),
+ MBCartVect(ray_vect),
+ t_enter, t_exit );
+}
+
static MBErrorCode intersect_children_with_elems(
MBAdaptiveKDTree* tool,
const MBRange& elems,
Modified: MOAB/trunk/MBAdaptiveKDTree.hpp
===================================================================
--- MOAB/trunk/MBAdaptiveKDTree.hpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/MBAdaptiveKDTree.hpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -393,6 +393,25 @@
//! if current node is root or back() will move to the
//! immediate sibling.
bool sibling_is_forward( ) const;
+
+ //! Find range of overlap between ray and leaf.
+ //!
+ //!\param ray_point Coordinates of start point of ray
+ //!\param ray_vect Directionion vector for ray such that
+ //! the ray is defined by r(t) = ray_point + t * ray_vect
+ //! for t > 0.
+ //!\param t_enter Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray entered the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\param t_exit Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray exited the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\return true if ray intersects leaf, false otherwise.
+ bool intersect_ray( const double ray_point[3],
+ const double ray_vect[3],
+ double& t_enter, double& t_exit ) const;
};
#endif
Modified: MOAB/trunk/MBBSPTree.cpp
===================================================================
--- MOAB/trunk/MBBSPTree.cpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/MBBSPTree.cpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -22,6 +22,7 @@
#include "MBGeomUtil.hpp"
#include "MBRange.hpp"
#include "MBInternals.hpp"
+#include "BSPTreePoly.hpp"
#include <assert.h>
#include <algorithm>
@@ -188,6 +189,12 @@
return moab()->tag_get_data( rootTag, &root_handle, 1, corners );
}
+MBErrorCode MBBSPTree::get_tree_box( MBEntityHandle root_handle,
+ double corners[24] )
+{
+ return moab()->tag_get_data( rootTag, &root_handle, 1, corners );
+}
+
MBErrorCode MBBSPTree::create_tree( MBEntityHandle& root_handle )
{
const double min[3] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL };
@@ -507,6 +514,13 @@
return tool()->get_split_plane( parent, plane );
}
+double MBBSPTreeIter::volume() const
+{
+ BSPTreePoly polyhedron;
+ MBErrorCode rval = calculate_polyhedron( polyhedron );
+ return MB_SUCCESS == rval ? polyhedron.volume() : -1.0;
+}
+
bool MBBSPTreeIter::is_sibling( const MBBSPTreeIter& other_leaf ) const
{
const size_t s = mStack.size();
@@ -543,6 +557,46 @@
return childVect[0] == handle();
}
+MBErrorCode MBBSPTreeIter::calculate_polyhedron( BSPTreePoly& poly_out ) const
+{
+ MBErrorCode rval;
+
+ assert( sizeof(MBCartVect) == 3*sizeof(double) );
+ MBCartVect corners[8];
+ rval = treeTool->get_tree_box( mStack.front(), corners[0].array() );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ rval = poly_out.set( corners );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ MBBSPTree::Plane plane;
+ std::vector<MBEntityHandle>::const_iterator i = mStack.begin();
+ std::vector<MBEntityHandle>::const_iterator here = mStack.end() - 1;
+ while (i != here) {
+ rval = treeTool->get_split_plane( *i, plane );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ childVect.clear();
+ rval = treeTool->moab()->get_child_meshsets( *i, childVect );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (childVect.size() != 2)
+ return MB_FAILURE;
+
+ ++i;
+ if (childVect[1] == *i)
+ plane.flip();
+
+ MBCartVect norm( plane.norm );
+ poly_out.cut_polyhedron( norm, plane.coeff );
+ }
+
+ return MB_SUCCESS;
+}
+
MBErrorCode MBBSPTreeBoxIter::initialize( MBBSPTree* tool_ptr,
MBEntityHandle root,
const double* point )
@@ -1154,6 +1208,11 @@
return MB_SUCCESS;
}
+MBErrorCode MBBSPTreeBoxIter::calculate_polyhedron( BSPTreePoly& poly_out ) const
+{
+ return poly_out.set( reinterpret_cast<const MBCartVect*>(leafCoords) );
+}
+
MBErrorCode MBBSPTree::leaf_containing_point( MBEntityHandle tree_root,
const double point[3],
MBEntityHandle& leaf_out )
@@ -1205,3 +1264,187 @@
return rval;
}
}
+
+
+
+template <typename PlaneIter> static inline
+bool ray_intersect_halfspaces( const MBCartVect& ray_pt,
+ const MBCartVect& ray_dir,
+ const PlaneIter& begin,
+ const PlaneIter& end,
+ double& t_enter,
+ double& t_exit )
+{
+ const double epsilon = 1e-12;
+
+ // begin with inifinite ray
+ t_enter = 0.0;
+ t_exit = std::numeric_limits<double>::infinity();
+
+ // cut ray such that we keep only the portion contained
+ // in each halfspace
+ for (PlaneIter i = begin; i != end; ++i) {
+ MBCartVect norm( i->norm );
+ double coeff = i->coeff;
+ double den = norm % ray_dir;
+ if (fabs(den) < epsilon) { // ray is parallel to plane
+ if (i->above( ray_pt.array() ))
+ return false; // ray entirely outside half-space
+ }
+ else {
+ double t_xsect = (-coeff - (norm % ray_pt)) / den;
+ // keep portion of ray/segment below plane
+ if (den > 0) {
+ if (t_xsect < t_exit)
+ t_exit = t_xsect;
+ }
+ else {
+ if (t_xsect > t_enter)
+ t_enter = t_xsect;
+ }
+ }
+ }
+
+ return t_exit >= t_enter;
+}
+
+class BoxPlaneIter {
+ int faceNum;
+ MBBSPTree::Plane facePlanes[6];
+
+ public:
+ BoxPlaneIter( const double coords[8][3] );
+ BoxPlaneIter( ) : faceNum(6) {} // initialize to 'end'
+ const MBBSPTree::Plane* operator->() const
+ { return facePlanes + faceNum; }
+ bool operator==( const BoxPlaneIter& other ) const
+ { return faceNum == other.faceNum; }
+ bool operator!=( const BoxPlaneIter& other ) const
+ { return faceNum != other.faceNum; }
+ BoxPlaneIter& operator++()
+ { ++faceNum; return *this; }
+};
+
+static const int box_face_corners[6][4] = { { 0, 1, 5, 4 },
+ { 1, 2, 6, 5 },
+ { 2, 3, 7, 6 },
+ { 3, 0, 4, 7 },
+ { 3, 2, 1, 0 },
+ { 4, 5, 6, 7 } };
+
+BoxPlaneIter::BoxPlaneIter( const double coords[8][3] )
+ : faceNum(0)
+{
+ // NOTE: In the case of a BSP tree, all sides of the
+ // leaf will planar.
+ assert( sizeof(MBCartVect) == sizeof(coords[0]) );
+ const MBCartVect* corners = reinterpret_cast<const MBCartVect*>(coords);
+ for (int i = 0; i < 6; ++i) {
+ const int* indices = box_face_corners[i];
+ MBCartVect v1 = corners[indices[1]] - corners[indices[0]];
+ MBCartVect v2 = corners[indices[3]] - corners[indices[0]];
+ MBCartVect n = v1 * v2;
+ facePlanes[i] = MBBSPTree::Plane( n.array(), -(n % corners[indices[2]]) );
+ }
+}
+
+
+bool MBBSPTreeBoxIter::intersect_ray( const double ray_point[3],
+ const double ray_vect[3],
+ double& t_enter, double& t_exit ) const
+{
+ BoxPlaneIter iter( this->leafCoords ), end;
+ return ray_intersect_halfspaces( MBCartVect(ray_point),
+ MBCartVect(ray_vect),
+ iter, end,
+ t_enter, t_exit );
+}
+
+class BSPTreePlaneIter {
+ MBBSPTree* toolPtr;
+ const MBEntityHandle* const pathToRoot;
+ int pathPos;
+ MBBSPTree::Plane tmpPlane;
+ std::vector<MBEntityHandle> tmpChildren;
+ public:
+ BSPTreePlaneIter( MBBSPTree* tool, const MBEntityHandle* path, int path_len )
+ : toolPtr(tool), pathToRoot(path), pathPos(path_len-1)
+ { operator++(); }
+ BSPTreePlaneIter() // initialize to 'end'
+ : toolPtr(0), pathToRoot(0), pathPos(-1) {}
+
+ const MBBSPTree::Plane* operator->() const
+ { return &tmpPlane; }
+ bool operator==( const BSPTreePlaneIter& other ) const
+ { return pathPos == other.pathPos; }
+ bool operator!=( const BSPTreePlaneIter& other ) const
+ { return pathPos != other.pathPos; }
+ BSPTreePlaneIter& operator++();
+};
+
+BSPTreePlaneIter& BSPTreePlaneIter::operator++()
+{
+ if (--pathPos < 0)
+ return *this;
+
+ MBEntityHandle prev = pathToRoot[pathPos+1];
+ MBEntityHandle curr = pathToRoot[pathPos];
+
+ MBErrorCode rval = toolPtr->get_split_plane( curr, tmpPlane );
+ if (MB_SUCCESS != rval) {
+ assert(false);
+ pathPos = 0;
+ return *this;
+ }
+
+ tmpChildren.clear();
+ rval = toolPtr->moab()->get_child_meshsets( curr, tmpChildren );
+ if (MB_SUCCESS != rval || tmpChildren.size() != 2) {
+ assert(false);
+ pathPos = 0;
+ return *this;
+ }
+
+ if (tmpChildren[1] == prev)
+ tmpPlane.flip();
+ return *this;
+}
+
+
+bool MBBSPTreeIter::intersect_ray( const double ray_point[3],
+ const double ray_vect[3],
+ double& t_enter, double& t_exit ) const
+{
+ // intersect with half-spaces defining tree
+ BSPTreePlaneIter iter1( tool(), &mStack[0], mStack.size() ), end1;
+ if (!ray_intersect_halfspaces( MBCartVect(ray_point),
+ MBCartVect(ray_vect),
+ iter1, end1,
+ t_enter, t_exit ))
+ return false;
+
+
+ // itersect with box bounding entire tree
+ double corners[8][3];
+ MBErrorCode rval = tool()->get_tree_box( mStack.front(), corners );
+ if (MB_SUCCESS != rval) {
+ assert(false);
+ return false;
+ }
+
+ BoxPlaneIter iter2( corners ), end2;
+ double t2_enter, t2_exit;
+ if (!ray_intersect_halfspaces( MBCartVect(ray_point),
+ MBCartVect(ray_vect),
+ iter2, end2,
+ t2_enter, t2_exit ))
+ return false;
+
+ // if intersect both box and halfspaces, check that
+ // two intersections overlap
+ if (t_enter < t2_enter)
+ t_enter = t2_enter;
+ if (t_exit > t2_exit)
+ t_exit = t2_exit;
+ return t_enter <= t_exit;
+}
Modified: MOAB/trunk/MBBSPTree.hpp
===================================================================
--- MOAB/trunk/MBBSPTree.hpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/MBBSPTree.hpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -31,6 +31,7 @@
class MBBSPTreeBoxIter;
class MBBSPTreeIter;
class MBRange;
+class BSPTreePoly;
class MBBSPTree
{
@@ -134,6 +135,10 @@
MBErrorCode get_tree_box( MBEntityHandle root_node,
double corner_coords[8][3] );
+ //! Get bounding box for entire tree
+ MBErrorCode get_tree_box( MBEntityHandle root_node,
+ double corner_coords[24] );
+
//! Set bounding box for entire tree
MBErrorCode set_tree_box( MBEntityHandle root_node,
const double box_min[3],
@@ -281,6 +286,28 @@
//! Get split plane that separates this node from its immediate sibling.
MBErrorCode get_parent_split_plane( MBBSPTree::Plane& plane ) const;
+ //! Get volume of leaf polyhedron
+ virtual double volume() const;
+
+ //! Find range of overlap between ray and leaf.
+ //!
+ //!\param ray_point Coordinates of start point of ray
+ //!\param ray_vect Directionion vector for ray such that
+ //! the ray is defined by r(t) = ray_point + t * ray_vect
+ //! for t > 0.
+ //!\param t_enter Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray entered the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\param t_exit Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray exited the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\return true if ray intersects leaf, false otherwise.
+ virtual bool intersect_ray( const double ray_point[3],
+ const double ray_vect[3],
+ double& t_enter, double& t_exit ) const;
+
//! Return true if thos node and the passed node share the
//! same immediate parent.
bool is_sibling( const MBBSPTreeIter& other_leaf ) const;
@@ -294,6 +321,9 @@
//! if current node is root or back() will move to the
//! immediate sibling.
bool sibling_is_forward( ) const;
+
+ //! Calculate the convex polyhedron bounding this leaf.
+ virtual MBErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
};
class MBBSPTreeBoxIter : public MBBSPTreeIter
@@ -370,6 +400,25 @@
//! test if a plane intersects the leaf box
bool intersects( const MBBSPTree::Plane& plane ) const;
+ //! Find range of overlap between ray and leaf.
+ //!
+ //!\param ray_point Coordinates of start point of ray
+ //!\param ray_vect Directionion vector for ray such that
+ //! the ray is defined by r(t) = ray_point + t * ray_vect
+ //! for t > 0.
+ //!\param t_enter Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray entered the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\param t_exit Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray exited the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\return true if ray intersects leaf, false otherwise.
+ bool intersect_ray( const double ray_point[3],
+ const double ray_vect[3],
+ double& t_enter, double& t_exit ) const;
+
//! Return the side of the box bounding this tree node
//! that is shared with the immediately adjacent sibling
//! (the tree node that shares a common parent node with
@@ -396,6 +445,9 @@
MBErrorCode get_neighbors( SideBits side,
std::vector<MBBSPTreeBoxIter>& results,
double epsilon = 0.0 ) const;
+
+ //! Calculate the convex polyhedron bounding this leaf.
+ MBErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
};
#endif
Modified: MOAB/trunk/MBGeomUtil.cpp
===================================================================
--- MOAB/trunk/MBGeomUtil.cpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/MBGeomUtil.cpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -156,6 +156,59 @@
return true;
}
+bool ray_box_intersect( const MBCartVect& box_min,
+ const MBCartVect& box_max,
+ const MBCartVect& ray_pt,
+ const MBCartVect& ray_dir,
+ double& t_enter, double& t_exit )
+{
+ const double epsilon = 1e-12;
+ double t1, t2;
+
+ // Use 'slabs' method from 13.6.1 of Akenine-Moller
+ t_enter = 0.0;
+ t_exit = std::numeric_limits<double>::infinity();
+
+ // Intersect with each pair of axis-aligned planes bounding
+ // opposite faces of the leaf box
+ bool ray_is_valid = false; // is ray direction vector zero?
+ for (int axis = 0; axis < 3; ++axis) {
+ if (fabs(ray_dir[axis]) < epsilon) { // ray parallel to planes
+ if (ray_pt[axis] >= box_min[axis] &&
+ ray_pt[axis] <= box_max[axis])
+ continue;
+ else
+ return false;
+ }
+
+ // find t values at which ray intersects each plane
+ ray_is_valid = true;
+ t1 = (box_min[axis] - ray_pt[axis]) / ray_dir[axis];
+ t2 = (box_max[axis] - ray_pt[axis]) / ray_dir[axis];
+
+ // t_enter = max( t_enter_x, t_enter_y, t_enter_z )
+ // t_exit = min( t_exit_x, t_exit_y, t_exit_z )
+ // where
+ // t_enter_x = min( t1_x, t2_x );
+ // t_exit_x = max( t1_x, t2_x )
+ if (t1 < t2) {
+ if (t_enter < t1)
+ t_enter = t1;
+ if (t_exit > t2)
+ t_exit = t2;
+ }
+ else {
+ if (t_enter < t2)
+ t_enter = t2;
+ if (t_exit > t1)
+ t_exit = t1;
+ }
+ }
+
+ return ray_is_valid && (t_enter <= t_exit);
+}
+
+
bool box_plane_overlap( const MBCartVect& normal,
double d,
MBCartVect min,
Modified: MOAB/trunk/MBGeomUtil.hpp
===================================================================
--- MOAB/trunk/MBGeomUtil.hpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/MBGeomUtil.hpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -75,6 +75,30 @@
double& t_out,
const double* ray_length = 0 );
+
+ //! Find range of overlap between ray and axis-aligned box.
+ //!
+ //!\param box_min Box corner with minimum coordinate values
+ //!\param box_max Box corner with minimum coordinate values
+ //!\param ray_pt Coordinates of start point of ray
+ //!\param ray_dir Directionion vector for ray such that
+ //! the ray is defined by r(t) = ray_point + t * ray_vect
+ //! for t > 0.
+ //!\param t_enter Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray entered the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\param t_exit Output: if return value is true, this value
+ //! is the parameter location along the ray at which
+ //! the ray exited the leaf. If return value is false,
+ //! then this value is undefined.
+ //!\return true if ray intersects leaf, false otherwise.
+bool ray_box_intersect( const MBCartVect& box_min,
+ const MBCartVect& box_max,
+ const MBCartVect& ray_pt,
+ const MBCartVect& ray_dir,
+ double& t_enter, double& t_exit );
+
/**\brief Test if plane intersects axis-aligned box
*
* Test for intersection between an unbounded plane and
Modified: MOAB/trunk/Makefile.am
===================================================================
--- MOAB/trunk/Makefile.am 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/Makefile.am 2009-09-25 15:34:25 UTC (rev 3168)
@@ -39,7 +39,8 @@
mesh_set_test \
cub_file_test \
exodus_test \
- mbcn_test
+ mbcn_test \
+ bsp_tree_poly_test
# merge_test \ # input files no longer exist?
# test_tag_server \ # fails
@@ -146,6 +147,7 @@
MBMeshSet.hpp \
MBOrientedBox.cpp \
MBOrientedBoxTreeTool.cpp \
+ BSPTreePoly.cpp \
MBRange.cpp \
MBRangeSeqIntersectIter.cpp \
MBRangeSeqIntersectIter.hpp \
@@ -258,6 +260,7 @@
MBOrientedBox.hpp \
MBOrientedBoxTreeTool.hpp \
MBParallelConventions.h \
+ BSPTreePoly.hpp \
MBRange.hpp \
MBReadUtilIface.hpp \
MBReaderIface.hpp \
@@ -368,6 +371,8 @@
mbcn_test_CPPFLAGS = -DSRCDIR=$(srcdir) # define anything to work around build issue
mbcn_test_LDADD =
+bsp_tree_poly_test_SOURCES = bsp_tree_poly_test.cpp
+
# Other files to clean up (e.g. output from tests)
DISTCLEANFILES = MBCN_FCDefs.h
MOSTLYCLEANFILES = dumped_acis.sat tree.h5m
Modified: MOAB/trunk/adaptive_kd_tree_tests.cpp
===================================================================
--- MOAB/trunk/adaptive_kd_tree_tests.cpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/adaptive_kd_tree_tests.cpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -1216,6 +1216,77 @@
CHECK( !iter.intersects( z6 ) );
}
+#define CHECK_RAY_XSECTS( PT, DIR, T_IN, T_OUT ) do { \
+ CHECK(iter.intersect_ray( (PT), (DIR), t_in, t_out )); \
+ CHECK_REAL_EQUAL( (T_IN), t_in, 1e-6 ); \
+ CHECK_REAL_EQUAL( (T_OUT), t_out, 1e-6 ); \
+ } while(false)
+
+void test_leaf_intersects_ray()
+{
+ MBErrorCode rval;
+ MBCore moab;
+ MBAdaptiveKDTree tool( &moab );
+ double t_in, t_out;
+
+ MBEntityHandle root;
+ const double min[3] = { -5, -4, -1 };
+ const double max[3] = { 1, 2, 3 };
+ rval = tool.create_tree( min, max, root );
+ CHECK_ERR(rval);
+
+ MBAdaptiveKDTreeIter iter;
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // start with point inside box
+ const double pt1[] = { 0, 0, 0 };
+ const double dir1[] = { 0.1, 0.1, 0.1 };
+ CHECK_RAY_XSECTS( pt1, dir1, 0, 10 );
+ const double dir2[] = { 5, 5, 5 };
+ CHECK_RAY_XSECTS( pt1, dir2, 0, 0.2 );
+ const double pxdir[] = { 1, 0, 0 };
+ CHECK_RAY_XSECTS( pt1, pxdir, 0, 1 );
+ const double nxdir[] = { -1, 0, 0 };
+ CHECK_RAY_XSECTS( pt1, nxdir, 0, 5 );
+ const double pydir[] = { 0, 1, 0 };
+ CHECK_RAY_XSECTS( pt1, pydir, 0, 2 );
+ const double nydir[] = { 0, -1, 0 };
+ CHECK_RAY_XSECTS( pt1, nydir, 0, 4 );
+ const double pzdir[] = { 0, 0, 1 };
+ CHECK_RAY_XSECTS( pt1, pzdir, 0, 3 );
+ const double nzdir[] = { 0, 0, -1 };
+ CHECK_RAY_XSECTS( pt1, nzdir, 0, 1 );
+
+ // point below box
+ const double pt2[] = { 0, 0, -2 };
+ CHECK_RAY_XSECTS( pt2, dir1, 10, 10 );
+ CHECK_RAY_XSECTS( pt2, dir2, 0.2, 0.2 );
+ CHECK(!iter.intersect_ray( pt2, pxdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, nxdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, pydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, nydir, t_in, t_out ));
+ CHECK_RAY_XSECTS( pt2, pzdir, 1, 5 );
+ CHECK(!iter.intersect_ray( pt2, nzdir, t_in, t_out ));
+
+ // point right of box
+ const double pt3[] = { 3, 0, 0 };
+ CHECK(!iter.intersect_ray( pt3, dir1, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, dir2, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, pxdir, t_in, t_out ));
+ CHECK_RAY_XSECTS( pt3, nxdir, 2, 8 );
+ CHECK(!iter.intersect_ray( pt3, pydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, nydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, pzdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, nzdir, t_in, t_out ));
+
+ // a few more complex test cases
+ const double dira[] = { -3, 0, 3 };
+ CHECK_RAY_XSECTS( pt3, dira, 2./3., 1.0 );
+ const double dirb[] = { -2, 0, 3.1 };
+ CHECK(!iter.intersect_ray( pt3, dirb, t_in, t_out ));
+}
+
int main()
{
int error_count = 0;
@@ -1229,5 +1300,6 @@
error_count += RUN_TEST(test_leaf_volume);
error_count += RUN_TEST(test_leaf_sibling);
error_count += RUN_TEST(test_leaf_intersects_plane);
+ error_count += RUN_TEST(test_leaf_intersects_ray);
return error_count;
}
Modified: MOAB/trunk/bsp_tree_test.cpp
===================================================================
--- MOAB/trunk/bsp_tree_test.cpp 2009-09-25 15:29:09 UTC (rev 3167)
+++ MOAB/trunk/bsp_tree_test.cpp 2009-09-25 15:34:25 UTC (rev 3168)
@@ -1,6 +1,9 @@
#include "MBCore.hpp"
#include "TestUtil.hpp"
#include "MBBSPTree.hpp"
+#include "MBCartVect.hpp"
+#include "BSPTreePoly.hpp"
+#include "MBRange.hpp"
#include <algorithm>
void test_set_plane();
@@ -13,8 +16,14 @@
void test_merge_leaf();
void test_box_iter_neighbors();
void test_leaf_sibling();
-void test_leaf_volume();
+void test_leaf_volume( bool box );
+void test_leaf_volume_box() { test_leaf_volume(true); }
+void test_leaf_volume_gen() { test_leaf_volume(false); }
void test_leaf_splits_intersects();
+void test_leaf_intersects_ray_common( bool box );
+void test_box_leaf_intersects_ray() { test_leaf_intersects_ray_common(true); }
+void test_gen_leaf_intersects_ray() { test_leaf_intersects_ray_common(false); }
+void test_leaf_polyhedron();
int main()
{
@@ -30,8 +39,12 @@
failures += RUN_TEST( test_merge_leaf );
failures += RUN_TEST( test_box_iter_neighbors );
failures += RUN_TEST( test_leaf_sibling );
- failures += RUN_TEST( test_leaf_volume );
+ failures += RUN_TEST( test_leaf_volume_box );
+ failures += RUN_TEST( test_leaf_volume_gen );
failures += RUN_TEST( test_leaf_splits_intersects );
+ failures += RUN_TEST( test_box_leaf_intersects_ray );
+ failures += RUN_TEST( test_gen_leaf_intersects_ray );
+ failures += RUN_TEST( test_leaf_polyhedron );
return failures;
}
@@ -1566,13 +1579,15 @@
CHECK( !iter1.sibling_is_forward() );
}
-void test_leaf_volume()
+void test_leaf_volume( bool box )
{
MBCore moab;
MBBSPTree tool( &moab );
MBErrorCode rval;
MBEntityHandle root;
- MBBSPTreeBoxIter iter;
+ MBBSPTreeBoxIter b_iter;
+ MBBSPTreeIter g_iter;
+ MBBSPTreeIter& iter = box ? b_iter : g_iter;
/* Build Tree
@@ -1701,3 +1716,386 @@
CHECK( iter.intersects( p ) );
}
+#define CHECK_RAY_XSECTS( PT, DIR, T_IN, T_OUT ) do { \
+ CHECK(iter.intersect_ray( (PT), (DIR), t_in, t_out )); \
+ CHECK_REAL_EQUAL( (T_IN), t_in, 1e-6 ); \
+ CHECK_REAL_EQUAL( (T_OUT), t_out, 1e-6 ); \
+ } while(false)
+
+void test_leaf_intersects_ray_common( bool box )
+{
+ MBErrorCode rval;
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ double t_in, t_out;
+
+ /** Start with only root box for initial testing **/
+
+ /* (1,5,-3) (2,5,-3)
+ o----o
+ /: / \
+ / : / \
+ (1,5,-1)o----o \
+ | : \ \
+ | : Y \ \
+ | : ^ \ \
+ | : | \ \
+ | : +-->X \ \ (6,1,-3)
+ | o./..........\.......o
+ | . L \ /
+ |. Z \ /
+ o--------------------o
+ (1,1,-1) (6,1,-1)
+ */
+ MBEntityHandle root;
+ const double corners[][3] = { { 1, 1, -3 },
+ { 6, 1, -3 },
+ { 2, 5, -3 },
+ { 1, 5, -3 },
+ { 1, 1, -1 },
+ { 6, 1, -1 },
+ { 2, 5, -1 },
+ { 1, 5, -1 } };
+ rval = tool.create_tree( corners, root );
+ CHECK_ERR(rval);
+
+ MBBSPTreeIter gen_iter;
+ MBBSPTreeBoxIter box_iter;
+ MBBSPTreeIter& iter = box ? static_cast<MBBSPTreeIter&>(box_iter) : gen_iter;
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // start with point inside box
+ const double pt1[] = { 3.5, 3, -2 };
+ const double dir1[] = { 0.1, 0.1, 0.1 };
+ CHECK_RAY_XSECTS( pt1, dir1, 0, 2.5 );
+ const double dir2[] = { 5, 5, 5 };
+ CHECK_RAY_XSECTS( pt1, dir2, 0, 0.05 );
+ const double pxdir[] = { 1, 0, 0 };
+ CHECK_RAY_XSECTS( pt1, pxdir, 0, 0.5 );
+ const double nxdir[] = { -1, 0, 0 };
+ CHECK_RAY_XSECTS( pt1, nxdir, 0, 2.5 );
+ const double pydir[] = { 0, 1, 0 };
+ CHECK_RAY_XSECTS( pt1, pydir, 0, 0.5 );
+ const double nydir[] = { 0, -1, 0 };
+ CHECK_RAY_XSECTS( pt1, nydir, 0, 2 );
+ const double pzdir[] = { 0, 0, 1 };
+ CHECK_RAY_XSECTS( pt1, pzdir, 0, 1 );
+ const double nzdir[] = { 0, 0, -1 };
+ CHECK_RAY_XSECTS( pt1, nzdir, 0, 1 );
+
+ // point below box
+ const double pt2[] = { 3.5, 3, -4 };
+ CHECK(!iter.intersect_ray( pt2, dir1, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, dir2, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, pxdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, nxdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, pydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt2, nydir, t_in, t_out ));
+ CHECK_RAY_XSECTS( pt2, pzdir, 1, 3 );
+ CHECK(!iter.intersect_ray( pt2, nzdir, t_in, t_out ));
+
+ // point right of box
+ const double pt3[] = { 7, 3, -2 };
+ CHECK(!iter.intersect_ray( pt3, dir1, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, dir2, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, pxdir, t_in, t_out ));
+ CHECK_RAY_XSECTS( pt3, nxdir, 3, 6 );
+ CHECK(!iter.intersect_ray( pt3, pydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, nydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, pzdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt3, nzdir, t_in, t_out ));
+
+ // a few more complex test cases
+ const double dira[] = { -2, -2, 0 };
+ CHECK_RAY_XSECTS( pt3, dira, 0.75, 1.0 );
+ const double dirb[] = { -1, -2.1, 0 };
+ CHECK(!iter.intersect_ray( pt3, dirb, t_in, t_out ));
+
+
+ /** Now split twice and test the bottom right corne **/
+
+ MBBSPTree::Plane Y3( MBBSPTree::Y, 3.0 );
+ rval = tool.split_leaf( iter, Y3 );
+ CHECK_ERR(rval);
+ MBBSPTree::Plane X2( MBBSPTree::X, 2.0 );
+ rval = tool.split_leaf( iter, X2 );
+ CHECK_ERR(rval);
+ rval = iter.step();
+ CHECK_ERR(rval);
+
+
+ /*
+ (2,3,-3)
+ o--------o (4,3,-3)
+ /: / \
+ / : / \
+ (2,3,-1)o--------o(4,3,-1)\ (6,1,-3)
+ | o.......\.......o
+ | . \ /
+ |. \ /
+ o---------------o
+ (2,1,-1) (6,1,-1)
+ */
+
+
+ // start with point inside box
+ const double pt4[] = { 4, 2, -2 };
+ CHECK_RAY_XSECTS( pt4, dir1, 0, 5 );
+ CHECK_RAY_XSECTS( pt4, dir2, 0, 0.1 );
+ CHECK_RAY_XSECTS( pt4, pxdir, 0, 1 );
+ CHECK_RAY_XSECTS( pt4, nxdir, 0, 2 );
+ CHECK_RAY_XSECTS( pt4, pydir, 0, 1 );
+ CHECK_RAY_XSECTS( pt4, nydir, 0, 1 );
+ CHECK_RAY_XSECTS( pt4, pzdir, 0, 1 );
+ CHECK_RAY_XSECTS( pt4, nzdir, 0, 1 );
+
+ // point below box
+ const double pt5[] = { 4, 2, -4 };
+ CHECK(!iter.intersect_ray( pt5, dir1, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt5, dir2, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt5, pxdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt5, nxdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt5, pydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt5, nydir, t_in, t_out ));
+ CHECK_RAY_XSECTS( pt5, pzdir, 1, 3 );
+ CHECK(!iter.intersect_ray( pt5, nzdir, t_in, t_out ));
+
+ // point right of box
+ const double pt6[] = { 7, 2, -2 };
+ CHECK(!iter.intersect_ray( pt6, dir1, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt6, dir2, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt6, pxdir, t_in, t_out ));
+ CHECK_RAY_XSECTS( pt6, nxdir, 2, 5 );
+ CHECK(!iter.intersect_ray( pt6, pydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt6, nydir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt6, pzdir, t_in, t_out ));
+ CHECK(!iter.intersect_ray( pt6, nzdir, t_in, t_out ));
+
+ // a few more complex test cases
+ const double dird[] = { -2, -2, 0 };
+ CHECK_RAY_XSECTS( pt6, dird, 0.5, 0.5 );
+ const double dire[] = { -3, -2, 0 };
+ CHECK_RAY_XSECTS( pt6, dire, 0.4, 0.5 );
+ const double dirf[] = { -2, -2.1, 0 };
+ CHECK(!iter.intersect_ray( pt6, dirf, t_in, t_out ));
+}
+
+static void box( const double pts[], int num_pts, double minpt[3], double maxpt[3] )
+{
+ minpt[0] = maxpt[0] = pts[0];
+ minpt[1] = maxpt[1] = pts[1];
+ minpt[2] = maxpt[2] = pts[2];
+ for (int i = 1; i < num_pts; ++i) {
+ for (int d = 0; d < 3; ++d) {
+ if (pts[3*i+d] < minpt[d])
+ minpt[d] = pts[3*i+d];
+ if (pts[3*i+d] > maxpt[d])
+ maxpt[d] = pts[3*i+d];
+ }
+ }
+}
+
+static MBEntityHandle build_tree( const double points[], int num_points, MBBSPTree& tool )
+{
+ // create points
+ MBErrorCode rval;
+ std::vector<MBEntityHandle> pts(num_points);
+ for (int i = 0; i < num_points; ++i) {
+ rval = tool.moab()->create_vertex( points + 3*i, pts[i] );
+ CHECK_ERR(rval);
+ }
+
+ // calculate bounding box of tree
+ double minpt[3], maxpt[3];
+ box( points, num_points, minpt, maxpt );
+
+ // create initial (1-node) tree
+ MBEntityHandle root;
+ rval = tool.create_tree( minpt, maxpt, root );
+ CHECK_ERR(rval);
+
+ MBBSPTreeIter iter;
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ rval = tool.moab()->add_entities( root, &pts[0], pts.size() );
+ CHECK_ERR(rval);
+
+ // build tree
+ std::vector<MBEntityHandle> left_pts, right_pts;
+ std::vector<MBCartVect> coords(num_points), tmp_coords;
+ std::vector<double> coeffs;
+ for (; MB_SUCCESS == rval; rval = iter.step()) {
+ pts.clear();
+ rval = tool.moab()->get_entities_by_handle( iter.handle(), pts );
+ CHECK_ERR(rval);
+ while (pts.size() > 1) {
+
+ coords.resize(pts.size());
+ rval = tool.moab()->get_coords( &pts[0], pts.size(), coords[0].array() );
+ CHECK_ERR(rval);
+
+ // find two points far apart apart
+ std::vector<MBCartVect>* ptr;
+ if (coords.size() < 10)
+ ptr = &coords;
+ else {
+ tmp_coords.resize(16);
+ MBCartVect pn, px;
+ box( coords[0].array(), coords.size(), pn.array(), px.array() );
+ tmp_coords[8] = pn;
+ tmp_coords[9] = MBCartVect( px[0], pn[1], pn[2] );
+ tmp_coords[10] = MBCartVect( px[0], px[1], pn[2] );
+ tmp_coords[11] = MBCartVect( pn[0], px[1], pn[2] );
+ tmp_coords[12] = MBCartVect( pn[0], pn[1], px[2] );
+ tmp_coords[13] = MBCartVect( px[0], pn[1], px[2] );
+ tmp_coords[14] = px;
+ tmp_coords[15] = MBCartVect( pn[0], px[1], px[2] );
+ for (int i = 0; i < 8; ++i) {
+ tmp_coords[i] = coords[0];
+ for (size_t j = 1; j < coords.size(); ++j)
+ if ((coords[j]-tmp_coords[i+8]).length_squared() <
+ (tmp_coords[i] - tmp_coords[i+8]).length_squared())
+ tmp_coords[i] = coords[j];
+ }
+ tmp_coords.resize(8);
+ ptr = &tmp_coords;
+ }
+
+ size_t pt1, pt2;
+ double lsqr = -1;
+ for (size_t i = 0; i < ptr->size(); ++i) {
+ for (size_t j = 0; j < ptr->size();++j) {
+ double ls = ((*ptr)[i] - (*ptr)[j]).length_squared();
+ if (ls > lsqr) {
+ lsqr = ls;
+ pt1 = i;
+ pt2 = j;
+ }
+ }
+ }
+
+ // if all points are coincident
+ if (lsqr <= 1e-12)
+ break;
+
+ // define normal orthogonal to line through two points
+ MBCartVect norm = (*ptr)[pt1] - (*ptr)[pt2];
+ norm.normalize();
+
+ // find mean position for plane
+ double coeff = 0.0;
+ for (size_t i = 0; i < coords.size(); ++i)
+ coeff -= norm % coords[i];
+ coeff /= coords.size();
+
+ // left/right sort points
+ left_pts.clear();
+ right_pts.clear();
+ for (size_t i = 0; i < coords.size(); ++i) {
+ double d = -(norm % coords[i]);
+ if (d >= coeff)
+ left_pts.push_back( pts[i] );
+ else
+ right_pts.push_back( pts[i] );
+ }
+
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( norm.array(), coeff ), left_pts, right_pts );
+ CHECK_ERR(rval);
+ CHECK( !left_pts.empty() && !right_pts.empty() );
+ pts.swap( left_pts );
+ }
+
+// printf("Leaf %d contains %d vertices: ", (int)ID_FROM_HANDLE(iter.handle()),
+// (int)(pts.size()) );
+// for (size_t i = 0; i < pts.size(); ++i)
+// printf( "%d, ", (int)ID_FROM_HANDLE(pts[i]));
+// printf("\n");
+ }
+
+ CHECK(rval == MB_ENTITY_NOT_FOUND);
+
+
+ // verify that tree is constructed correctly
+ for (int i = 0; i < num_points; ++i) {
+ MBCartVect pt( points + 3*i );
+ MBEntityHandle leaf;
+ rval = tool.leaf_containing_point( root, pt.array(), leaf );
+ CHECK_ERR(rval);
+ MBRange ents;
+ rval = tool.moab()->get_entities_by_handle( leaf, ents );
+ CHECK_ERR(rval);
+ bool found = false;
+ for (MBRange::iterator j = ents.begin(); j != ents.end(); ++j) {
+ MBCartVect ent_coords;
+ rval = tool.moab()->get_coords( &*j, 1, ent_coords.array() );
+ CHECK_ERR(rval);
+ if ((pt - ent_coords).length_squared() < 1e-6)
+ found = true;
+ }
+ CHECK(found);
+ }
+
+ return root;
+}
+
+
+void test_leaf_polyhedron()
+{
+ // array of 20 points used to construct tree
+ static const double points[] = {
+ 7, 6, 3,
+ 5, 3, 5,
+ 9, 2, 6,
+ 7, 2, 1,
+ 3, 9, 0,
+ 6, 0, 6,
+ 1, 6, 2,
+ 9, 7, 8,
+ 2, 0, 2,
+ 5, 7, 3,
+ 2, 2, 9,
+ 7, 9, 8,
+ 1, 6, 3,
+ 3, 9, 2,
+ 4, 9, 1,
+ 4, 8, 7,
+ 3, 0, 5,
+ 0, 1, 6,
+ 2, 3, 6,
+ 1, 6, 0};
+ const int num_pts = sizeof(points)/(3*sizeof(double));
+
+ MBErrorCode rval;
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBEntityHandle root = build_tree( points, num_pts, tool );
+
+ MBBSPTreeIter iter;
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ std::vector<MBEntityHandle> pts;
+ std::vector<MBCartVect> coords;
+
+ for (; rval == MB_SUCCESS; rval = iter.step()) {
+ BSPTreePoly poly;
+ rval = iter.calculate_polyhedron( poly );
+ CHECK_ERR(rval);
+
+ CHECK( poly.is_valid() );
+ CHECK( poly.volume() > 0.0 );
+
+ pts.clear();
+ rval = tool.moab()->get_entities_by_handle( iter.handle(), pts );
+ CHECK_ERR(rval);
+ CHECK( !pts.empty() );
+ coords.resize(pts.size());
+ rval = tool.moab()->get_coords( &pts[0], pts.size(), coords[0].array() );
+ CHECK_ERR(rval);
+
+ for (size_t i = 0; i < pts.size(); ++i)
+ CHECK( poly.is_point_contained( coords[i] ) );
+ }
+}
More information about the moab-dev
mailing list