[MOAB-dev] r1712 - MOAB/trunk
kraftche at mcs.anl.gov
kraftche at mcs.anl.gov
Thu Mar 27 17:49:54 CDT 2008
Author: kraftche
Date: 2008-03-27 17:49:54 -0500 (Thu, 27 Mar 2008)
New Revision: 1712
Modified:
MOAB/trunk/MBAdaptiveKDTree.cpp
MOAB/trunk/MBAdaptiveKDTree.hpp
Log:
fix potentially very slow closest_triangle query
Modified: MOAB/trunk/MBAdaptiveKDTree.cpp
===================================================================
--- MOAB/trunk/MBAdaptiveKDTree.cpp 2008-03-27 18:55:08 UTC (rev 1711)
+++ MOAB/trunk/MBAdaptiveKDTree.cpp 2008-03-27 22:49:54 UTC (rev 1712)
@@ -1372,6 +1372,118 @@
return MB_SUCCESS;
}
+static MBErrorCode closest_to_triangles( MBInterface* moab,
+ const MBRange& tris,
+ const MBCartVect& from,
+ double& shortest_dist_sqr,
+ MBCartVect& closest_pt,
+ MBEntityHandle& closest_tri )
+{
+ MBErrorCode rval;
+ MBCartVect pos, diff, verts[3];
+ const MBEntityHandle* conn;
+ int len;
+
+ for (MBRange::iterator i = tris.begin(); i != tris.end(); ++i) {
+ rval = moab->get_connectivity( *i, conn, len );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ rval = moab->get_coords( conn, 3, verts[0].array() );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ MBGeomUtil::closest_location_on_tri( from, verts, pos );
+ diff = pos - from;
+ double dist_sqr = diff % diff;
+ if (dist_sqr < shortest_dist_sqr) {
+ // new closest location
+ shortest_dist_sqr = dist_sqr;
+ closest_pt = pos;
+ closest_tri = *i;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+
+static MBErrorCode closest_to_triangles( MBInterface* moab,
+ MBEntityHandle set_handle,
+ const MBCartVect& from,
+ double& shortest_dist_sqr,
+ MBCartVect& closest_pt,
+ MBEntityHandle& closest_tri )
+{
+ MBErrorCode rval;
+ MBRange tris;
+
+ rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ return closest_to_triangles( moab, tris, from, shortest_dist_sqr, closest_pt, closest_tri );
+}
+
+MBErrorCode MBAdaptiveKDTree::find_close_triangle( MBEntityHandle root,
+ const double from[3],
+ double pt[3],
+ MBEntityHandle& triangle )
+{
+ MBErrorCode rval;
+ MBRange tris;
+ Plane split;
+ std::vector<MBEntityHandle> stack;
+ std::vector<MBEntityHandle> children(2);
+ stack.reserve(30);
+ stack.push_back( root );
+
+ while (!stack.empty()) {
+ MBEntityHandle node = stack.back();
+ stack.pop_back();
+
+ for (;;) { // loop until we find a leaf
+
+ children.clear();
+ rval = moab()->get_child_meshsets( node, children );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // loop termination criterion
+ if (children.empty())
+ break;
+
+ // if not a leaf, get split plane
+ rval = get_split_plane( node, split );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // continue down the side that contains the point,
+ // and push the other side onto the stack in case
+ // we need to check it later.
+ int rs = split.right_side( from );
+ node = children[rs];
+ stack.push_back( children[1-rs] );
+ }
+
+ // We should now be at a leaf.
+ // If it has some triangles, we're done.
+ // If not, continue searching for another leaf.
+ rval = moab()->get_entities_by_type( node, MBTRI, tris );
+ if (!tris.empty()) {
+ double dist_sqr = HUGE_VAL;
+ MBCartVect point(pt);
+ rval = closest_to_triangles( moab(), tris, MBCartVect(from), dist_sqr, point, triangle );
+ point.get(pt);
+ return rval;
+ }
+ }
+
+ // If we got here, then we traversed the entire tree
+ // and all the leaves were empty.
+ return MB_ENTITY_NOT_FOUND;
+}
+
/** Find the triangles in a set that are closer to the input
* position than any triangles in the 'closest_tris' list.
*
@@ -1479,46 +1591,6 @@
return MB_SUCCESS;
}
-static MBErrorCode closest_to_triangles( MBInterface* moab,
- MBEntityHandle set_handle,
- const MBCartVect& from,
- double& shortest_dist_sqr,
- MBCartVect& closest_pt,
- MBEntityHandle& closest_tri )
-{
- MBErrorCode rval;
- MBRange tris;
- MBCartVect pos, diff, verts[3];
- const MBEntityHandle* conn;
- int len;
-
- rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
- if (MB_SUCCESS != rval)
- return rval;
-
- for (MBRange::iterator i = tris.begin(); i != tris.end(); ++i) {
- rval = moab->get_connectivity( *i, conn, len );
- if (MB_SUCCESS != rval)
- return rval;
-
- rval = moab->get_coords( conn, 3, verts[0].array() );
- if (MB_SUCCESS != rval)
- return rval;
-
- MBGeomUtil::closest_location_on_tri( from, verts, pos );
- diff = pos - from;
- double dist_sqr = diff % diff;
- if (dist_sqr < shortest_dist_sqr) {
- // new closest location
- shortest_dist_sqr = dist_sqr;
- closest_pt = pos;
- closest_tri = *i;
- }
- }
-
- return MB_SUCCESS;
-}
-
MBErrorCode MBAdaptiveKDTree::closest_triangle( MBEntityHandle tree_root,
const double from_coords[3],
double closest_point_out[3],
@@ -1533,15 +1605,9 @@
// Find the leaf containing the input point
// This search does not take into account any bounding box for the
// tree, so it always returns one leaf.
- MBEntityHandle leaf;
- rval = leaf_containing_point( tree_root, from_coords, leaf );
+ rval = find_close_triangle( tree_root, from_coords, closest_pt.array(), triangle_out );
if (MB_SUCCESS != rval) return rval;
- // Find the closest triangle(s) in the leaf containing the point
- rval = closest_to_triangles( moab(), leaf, from, shortest_dist_sqr,
- closest_pt, triangle_out );
- if (MB_SUCCESS != rval) return rval;
-
// Find any other leaves for which the bounding box is within
// the same distance from the input point as the current closest
// point is.
Modified: MOAB/trunk/MBAdaptiveKDTree.hpp
===================================================================
--- MOAB/trunk/MBAdaptiveKDTree.hpp 2008-03-27 18:55:08 UTC (rev 1711)
+++ MOAB/trunk/MBAdaptiveKDTree.hpp 2008-03-27 22:49:54 UTC (rev 1712)
@@ -25,6 +25,7 @@
#include <string>
#include <vector>
+#include <math.h>
class MBAdaptiveKDTreeIter;
class MBInterface;
@@ -51,6 +52,19 @@
struct Plane {
double coord; //!< Location of plane as coordinate on normal axis
int norm; //!< The principal axis that is the normal of the plane;
+
+ /** return true if point is below/to the left of the split plane */
+ bool left_side( const double point[3] ) {
+ return point[norm] < coord;
+ }
+ /** return true if point is abve/to the right of the split plane */
+ bool right_side( const double point[3] ) {
+ return point[norm] > coord;
+ }
+ /** return distance from point to plane */
+ double distance( const double point[3] ) const {
+ return fabs(point[norm] - coord);
+ }
};
//! Get split plane for tree node
@@ -201,6 +215,13 @@
const double distance,
std::vector<MBEntityHandle>& leaves_out );
+private:
+
+ /**\brief find a triangle near the input point */
+ MBErrorCode find_close_triangle( MBEntityHandle root,
+ const double from_point[3],
+ double pt[3],
+ MBEntityHandle& triangle );
};
More information about the moab-dev
mailing list