Branch data Line data Source code
1 : : /**\file AdaptiveKDTree.hpp
2 : : * \class moab::AdaptiveKDTree
3 : : * \brief Adaptive KD tree, for sorting and searching entities spatially
4 : : */
5 : :
6 : : #ifndef MOAB_ADAPTIVE_KD_TREE_HPP
7 : : #define MOAB_ADAPTIVE_KD_TREE_HPP
8 : :
9 : : #include "moab/Types.hpp"
10 : : #include "moab/Tree.hpp"
11 : :
12 : : #include <string>
13 : : #include <vector>
14 : : #include <math.h>
15 : :
16 : : namespace moab
17 : : {
18 : :
19 : : class AdaptiveKDTreeIter;
20 : : class Interface;
21 : : class Range;
22 : :
23 : : class AdaptiveKDTree : public Tree
24 : : {
25 : : public:
26 : : AdaptiveKDTree( Interface* iface );
27 : :
28 : : /** \brief Constructor (build the tree on construction)
29 : : * Construct a tree object, and build the tree with entities input. See comments
30 : : * for build_tree() for detailed description of arguments.
31 : : * \param iface MOAB instance
32 : : * \param entities Entities to build tree around
33 : : * \param tree_root Root set for tree (see function description)
34 : : * \param opts Options for tree (see function description)
35 : : */
36 : : AdaptiveKDTree( Interface* iface, const Range& entities, EntityHandle* tree_root_set = NULL,
37 : : FileOptions* opts = NULL );
38 : :
39 : : ~AdaptiveKDTree();
40 : :
41 : : /** \brief Parse options for tree creation
42 : : * \param options Options passed in by application
43 : : * \return Failure is returned if any options were passed in and not interpreted; could mean
44 : : * inappropriate options for a particular tree type
45 : : */
46 : : ErrorCode parse_options( FileOptions& options );
47 : :
48 : : /** Build the tree
49 : : * Build a tree with the entities input. If a non-NULL tree_root_set pointer is input,
50 : : * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct
51 : : * a new root set and pass its handle back in *tree_root_set. Options vary by tree type;
52 : : * see Tree.hpp for common options; options specific to AdaptiveKDTree:
53 : : * SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
54 : : * PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
55 : : * for possible values; default = 1 (SUBDIVISION_SNAP)
56 : : * \param entities Entities with which to build the tree
57 : : * \param tree_root Root set for tree (see function description)
58 : : * \param opts Options for tree (see function description)
59 : : * \return Error is returned only on build failure
60 : : */
61 : : virtual ErrorCode build_tree( const Range& entities, EntityHandle* tree_root_set = NULL,
62 : : FileOptions* options = NULL );
63 : :
64 : : //! Reset the tree, optionally checking we have the right root
65 : : virtual ErrorCode reset_tree();
66 : :
67 : : /** \brief Get leaf containing input position.
68 : : *
69 : : * Does not take into account global bounding box of tree.
70 : : * - Therefore there is always one leaf containing the point.
71 : : * - If caller wants to account for global bounding box, then
72 : : * caller can test against that box and not call this method
73 : : * at all if the point is outside the box, as there is no leaf
74 : : * containing the point in that case.
75 : : * \param point Point to be located in tree
76 : : * \param leaf_out Leaf containing point
77 : : * \param iter_tol Tolerance for convergence of point search
78 : : * \param inside_tol Tolerance for inside element calculation
79 : : * \param multiple_leaves Some tree types can have multiple leaves containing a point;
80 : : * if non-NULL, this parameter is returned true if multiple leaves contain
81 : : * the input point
82 : : * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
83 : : * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
84 : : */
85 : : virtual ErrorCode point_search( const double* point, EntityHandle& leaf_out, const double iter_tol = 1.0e-10,
86 : : const double inside_tol = 1.0e-6, bool* multiple_leaves = NULL,
87 : : EntityHandle* start_node = NULL, CartVect* params = NULL );
88 : :
89 : : /** \brief Get leaf containing input position.
90 : : *
91 : : * Does not take into account global bounding box of tree.
92 : : * - Therefore there is always one leaf containing the point.
93 : : * - If caller wants to account for global bounding box, then
94 : : * caller can test against that box and not call this method
95 : : * at all if the point is outside the box, as there is no leaf
96 : : * containing the point in that case.
97 : : * \param point Point to be located in tree
98 : : * \param leaf_it Iterator to leaf containing point
99 : : * \param iter_tol Tolerance for convergence of point search
100 : : * \param inside_tol Tolerance for inside element calculation
101 : : * \param multiple_leaves Some tree types can have multiple leaves containing a point;
102 : : * if non-NULL, this parameter is returned true if multiple leaves contain
103 : : * the input point
104 : : * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
105 : : * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
106 : : */
107 : : ErrorCode point_search( const double* point, AdaptiveKDTreeIter& leaf_it, const double iter_tol = 1.0e-10,
108 : : const double inside_tol = 1.0e-6, bool* multiple_leaves = NULL,
109 : : EntityHandle* start_node = NULL );
110 : :
111 : : /** \brief Find all leaves within a given distance from point
112 : : * If dists_out input non-NULL, also returns distances from each leaf; if
113 : : * point i is inside leaf, 0 is given as dists_out[i].
114 : : * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities
115 : : * in tree nodes and return containing entities in leaves_out. In those cases, if params_out
116 : : * is also non-NULL, will return parameters in those elements in that vector.
117 : : * \param point Point to be located in tree
118 : : * \param distance Distance within which to query
119 : : * \param leaves_out Leaves within distance or containing point
120 : : * \param iter_tol Tolerance for convergence of point search
121 : : * \param inside_tol Tolerance for inside element calculation
122 : : * \param dists_out If non-NULL, will contain distsances to leaves
123 : : * \param params_out If non-NULL, will contain parameters of the point in the ents in leaves_out
124 : : * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
125 : : */
126 : : virtual ErrorCode distance_search( const double* point, const double distance,
127 : : std::vector< EntityHandle >& leaves_out, const double iter_tol = 1.0e-10,
128 : : const double inside_tol = 1.0e-6, std::vector< double >* dists_out = NULL,
129 : : std::vector< CartVect >* params_out = NULL, EntityHandle* start_node = NULL );
130 : :
131 : : ErrorCode get_info( EntityHandle root, double min[3], double max[3], unsigned int& dep );
132 : :
133 : : //! Enumeriate split plane directions
134 : : enum Axis
135 : : {
136 : : X = 0,
137 : : Y = 1,
138 : : Z = 2
139 : : };
140 : :
141 : : //! Split plane
142 : : struct Plane
143 : : {
144 : : double coord; //!< Location of plane as coordinate on normal axis
145 : : int norm; //!< The principal axis that is the normal of the plane;
146 : :
147 : : /** return true if point is below/to the left of the split plane */
148 : : bool left_side( const double point[3] )
149 : : {
150 : : return point[norm] < coord;
151 : : }
152 : : /** return true if point is above/to the right of the split plane */
153 : 161 : bool right_side( const double point[3] )
154 : : {
155 : 161 : return point[norm] > coord;
156 : : }
157 : : /** return distance from point to plane */
158 : : double distance( const double point[3] ) const
159 : : {
160 : : return fabs( point[norm] - coord );
161 : : }
162 : : };
163 : :
164 : : //! Get split plane for tree node
165 : : ErrorCode get_split_plane( EntityHandle node, Plane& plane );
166 : :
167 : : //! Set split plane for tree node
168 : : ErrorCode set_split_plane( EntityHandle node, const Plane& plane );
169 : :
170 : : //! Get iterator for tree
171 : : ErrorCode get_tree_iterator( EntityHandle tree_root, AdaptiveKDTreeIter& result );
172 : :
173 : : //! Get iterator at right-most ('last') leaf.
174 : : ErrorCode get_last_iterator( EntityHandle tree_root, AdaptiveKDTreeIter& result );
175 : :
176 : : //! Get iterator for tree or subtree
177 : : ErrorCode get_sub_tree_iterator( EntityHandle tree_root, const double box_min[3], const double box_max[3],
178 : : AdaptiveKDTreeIter& result );
179 : :
180 : : //! Split leaf of tree
181 : : //! Updates iterator location to point to first new leaf node.
182 : : ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane );
183 : :
184 : : //! Split leaf of tree
185 : : //! Updates iterator location to point to first new leaf node.
186 : : ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, EntityHandle& left_child, EntityHandle& right_child );
187 : : //! Split leaf of tree
188 : : //! Updates iterator location to point to first new leaf node.
189 : : ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, const Range& left_entities,
190 : : const Range& right_entities );
191 : :
192 : : //! Split leaf of tree
193 : : //! Updates iterator location to point to first new leaf node.
194 : : ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, const std::vector< EntityHandle >& left_entities,
195 : : const std::vector< EntityHandle >& right_entities );
196 : :
197 : : //! Merge the leaf pointed to by the current iterator with it's
198 : : //! sibling. If the sibling is not a leaf, multiple merges may
199 : : //! be done.
200 : : ErrorCode merge_leaf( AdaptiveKDTreeIter& iter );
201 : :
202 : : //! Find triangle closest to input position.
203 : : //!\param from_coords The input position to test against
204 : : //!\param closest_point_out The closest point on the set of triangles in the tree
205 : : //!\param triangle_out The triangle closest to the input position
206 : : ErrorCode closest_triangle( EntityHandle tree_root, const double from_coords[3], double closest_point_out[3],
207 : : EntityHandle& triangle_out );
208 : :
209 : : ErrorCode sphere_intersect_triangles( EntityHandle tree_root, const double center[3], double radius,
210 : : std::vector< EntityHandle >& triangles );
211 : :
212 : : ErrorCode ray_intersect_triangles( EntityHandle tree_root, const double tolerance, const double ray_unit_dir[3],
213 : : const double ray_base_pt[3], std::vector< EntityHandle >& triangles_out,
214 : : std::vector< double >& distance_out, int result_count_limit = 0,
215 : : double distance_limit = -1.0 );
216 : :
217 : : ErrorCode compute_depth( EntityHandle root, unsigned int& min_depth, unsigned int& max_depth );
218 : :
219 : : //! methods for selecting candidate split planes
220 : : enum CandidatePlaneSet
221 : : {
222 : : //! Candidiate planes at evenly spaced intervals
223 : : SUBDIVISION = 0,
224 : : //! Like SUBDIVISION, except snap to closest vertex coordinate
225 : : SUBDIVISION_SNAP, // = 1
226 : : //! Median vertex coodinate values
227 : : VERTEX_MEDIAN, // = 2
228 : : //! Random sampling of vertex coordinate values
229 : : VERTEX_SAMPLE // = 3
230 : : };
231 : :
232 : : //! print various things about this tree
233 : : virtual ErrorCode print();
234 : :
235 : : private:
236 : : friend class AdaptiveKDTreeIter;
237 : :
238 : : ErrorCode init();
239 : :
240 : : /**\brief find a triangle near the input point */
241 : : ErrorCode find_close_triangle( EntityHandle root, const double from_point[3], double pt[3],
242 : : EntityHandle& triangle );
243 : :
244 : : ErrorCode make_tag( Interface* iface, std::string name, TagType storage, DataType type, int count,
245 : : void* default_val, Tag& tag_handle, std::vector< Tag >& created_tags );
246 : :
247 : : ErrorCode intersect_children_with_elems( const Range& elems, AdaptiveKDTree::Plane plane, double eps,
248 : : CartVect box_min, CartVect box_max, Range& left_tris, Range& right_tris,
249 : : Range& both_tris, double& metric_value );
250 : :
251 : : ErrorCode best_subdivision_snap_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
252 : : Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
253 : : std::vector< double >& tmp_data, double eps );
254 : :
255 : : ErrorCode best_subdivision_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
256 : : Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
257 : : double eps );
258 : :
259 : : ErrorCode best_vertex_median_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
260 : : Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
261 : : std::vector< double >& coords, double eps );
262 : :
263 : : ErrorCode best_vertex_sample_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
264 : : Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
265 : : std::vector< double >& coords, std::vector< EntityHandle >& indices,
266 : : double eps );
267 : :
268 : : static const char* treeName;
269 : :
270 : : Tag planeTag, axisTag;
271 : :
272 : : unsigned splitsPerDir;
273 : :
274 : : CandidatePlaneSet planeSet;
275 : :
276 : : bool spherical;
277 : : double radius;
278 : : };
279 : :
280 : : //! Iterate over leaves of an adaptive kD-tree
281 [ # # ][ + + ]: 205 : class AdaptiveKDTreeIter
[ + - ][ + + ]
282 : : {
283 : : public:
284 : : enum Direction
285 : : {
286 : : LEFT = 0,
287 : : RIGHT = 1
288 : : };
289 : :
290 : : private:
291 : : struct StackObj
292 : : {
293 : 5700 : StackObj( EntityHandle e, double c ) : entity( e ), coord( c ) {}
294 : 11404 : StackObj() : entity( 0 ), coord( 0.0 ) {}
295 : : EntityHandle entity; //!< handle for tree node
296 : : double coord; //!< box coordinate of parent
297 : : };
298 : :
299 : : enum
300 : : {
301 : : BMIN = 0,
302 : : BMAX = 1
303 : : }; //!< indices into mBox and child list
304 : :
305 : : CartVect mBox[2]; //!< min and max corners of bounding box
306 : : AdaptiveKDTree* treeTool; //!< tool for tree
307 : : std::vector< StackObj > mStack; //!< stack storing path through tree
308 : : mutable std::vector< EntityHandle > childVect; //!< temporary storage of child handles
309 : :
310 : : //! Descend tree to left most leaf from current position
311 : : //! No-op if at leaf.
312 : : ErrorCode step_to_first_leaf( Direction direction );
313 : :
314 : : friend class AdaptiveKDTree;
315 : :
316 : : public:
317 [ + + ][ + - ]: 153 : AdaptiveKDTreeIter() : treeTool( 0 ), childVect( 2 ) {}
318 : :
319 : : ErrorCode initialize( AdaptiveKDTree* tool, EntityHandle root, const double box_min[3], const double box_max[3],
320 : : Direction direction );
321 : :
322 : 10225 : AdaptiveKDTree* tool() const
323 : : {
324 : 10225 : return treeTool;
325 : : }
326 : :
327 : : //! Get handle for current leaf
328 : 28152 : EntityHandle handle() const
329 : : {
330 : 28152 : return mStack.back().entity;
331 : : }
332 : :
333 : : //! Get min corner of axis-aligned box for current leaf
334 : 4808 : const double* box_min() const
335 : : {
336 : 4808 : return mBox[BMIN].array();
337 : : }
338 : :
339 : : //! Get max corner of axis-aligned box for current leaf
340 : 4808 : const double* box_max() const
341 : : {
342 : 4808 : return mBox[BMAX].array();
343 : : }
344 : :
345 : 8 : double volume() const
346 : : {
347 : 8 : return ( mBox[BMAX][0] - mBox[BMIN][0] ) * ( mBox[BMAX][1] - mBox[BMIN][1] ) *
348 : 8 : ( mBox[BMAX][2] - mBox[BMIN][2] );
349 : : }
350 : :
351 : : //! test if a plane intersects the leaf box
352 : 18 : bool intersects( const AdaptiveKDTree::Plane& plane ) const
353 : : {
354 [ + + ][ + + ]: 18 : return mBox[BMIN][plane.norm] <= plane.coord && mBox[BMAX][plane.norm] >= plane.coord;
355 : : }
356 : :
357 : : //! Get depth in tree. root is at depth of 1.
358 : 6466 : unsigned depth() const
359 : : {
360 : 6466 : return mStack.size();
361 : : }
362 : :
363 : : //! Advance the iterator either left or right in the tree
364 : : //! Note: stepping past the end of the tree will invalidate
365 : : //! the iterator. It will *not* be work step the
366 : : //! other direction.
367 : : ErrorCode step( Direction direction );
368 : :
369 : : //! Advance to next leaf
370 : : //! Returns MB_ENTITY_NOT_FOUND if at end.
371 : : //! Note: steping past the end of the tree will invalidate
372 : : //! the iterator. Calling back() will not work.
373 : 5499 : ErrorCode step()
374 : : {
375 : 5499 : return step( RIGHT );
376 : : }
377 : :
378 : : //! Move back to previous leaf
379 : : //! Returns MB_ENTITY_NOT_FOUND if at beginning.
380 : : //! Note: steping past the start of the tree will invalidate
381 : : //! the iterator. Calling step() will not work.
382 : 191 : ErrorCode back()
383 : : {
384 : 191 : return step( LEFT );
385 : : }
386 : :
387 : : //! Return the side of the box bounding this tree node
388 : : //! that is shared with the immediately adjacent sibling
389 : : //! (the tree node that shares a common parent node with
390 : : //! this node in the binary tree.)
391 : : //!
392 : : //!\param axis_out The principal axis orthogonal to the side of the box
393 : : //!\param neg_out true if the side of the box is toward the decreasing
394 : : //! direction of the principal axis indicated by axis_out,
395 : : //! false if it is toward the increasing direction.
396 : : //!\return MB_ENTITY_NOT FOUND if root node.
397 : : //! MB_FAILURE if internal error.
398 : : //! MB_SUCCESS otherwise.
399 : : ErrorCode sibling_side( AdaptiveKDTree::Axis& axis_out, bool& neg_out ) const;
400 : :
401 : : //! Get adjacent leaf nodes on side indicated by norm and neg.
402 : : //!
403 : : //! E.g. if norm == X and neg == true, then get neighbor(s)
404 : : //! adjacent to the side of the box contained in the plane
405 : : //! with normal to the X axis and with the x coordinate equal
406 : : //! to the minimum x of the bounding box.
407 : : //!
408 : : //! E.g. if norm == Y and neg == false, then get neighbor(s)
409 : : //! adjacent to the side of the box with y = maximum y of bounding box.
410 : : //!
411 : : //!\param norm Normal vector for box side (X, Y, or Z)
412 : : //!\param neg Which of two planes with norm (true->smaller coord,
413 : : //! false->larger coord)
414 : : //!\param results List to which to append results. This function does
415 : : //! *not* clear existing values in list.
416 : : //!\param epsilon Tolerance on overlap. A positive value E will
417 : : //! result in nodes that are separated by as much as E
418 : : //! to be considered touching. A negative value -E will
419 : : //! cause leaves that do not overlap by at least E to be
420 : : //! considered non-overlapping. Amongst other things,
421 : : //! this value can be used to control whether or not
422 : : //! leaves adjacent at only their edges or corners are
423 : : //! returned.
424 : : ErrorCode get_neighbors( AdaptiveKDTree::Axis norm, bool neg, std::vector< AdaptiveKDTreeIter >& results,
425 : : double epsilon = 0.0 ) const;
426 : :
427 : : //! Get split plane that separates this node from its immediate sibling.
428 : : ErrorCode get_parent_split_plane( AdaptiveKDTree::Plane& plane ) const;
429 : :
430 : : //! Return true if thos node and the passed node share the
431 : : //! same immediate parent.
432 : : bool is_sibling( const AdaptiveKDTreeIter& other_leaf ) const;
433 : :
434 : : //! Return true if thos node and the passed node share the
435 : : //! same immediate parent.
436 : : bool is_sibling( EntityHandle other_leaf ) const;
437 : :
438 : : //! Returns true if calling step() will advance to the
439 : : //! immediate sibling of the current node. Returns false
440 : : //! if current node is root or back() will move to the
441 : : //! immediate sibling.
442 : : bool sibling_is_forward() const;
443 : :
444 : : //! Find range of overlap between ray and leaf.
445 : : //!
446 : : //!\param ray_point Coordinates of start point of ray
447 : : //!\param ray_vect Directionion vector for ray such that
448 : : //! the ray is defined by r(t) = ray_point + t * ray_vect
449 : : //! for t > 0.
450 : : //!\param t_enter Output: if return value is true, this value
451 : : //! is the parameter location along the ray at which
452 : : //! the ray entered the leaf. If return value is false,
453 : : //! then this value is undefined.
454 : : //!\param t_exit Output: if return value is true, this value
455 : : //! is the parameter location along the ray at which
456 : : //! the ray exited the leaf. If return value is false,
457 : : //! then this value is undefined.
458 : : //!\return true if ray intersects leaf, false otherwise.
459 : : bool intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter, double& t_exit ) const;
460 : : };
461 : :
462 : 29 : inline ErrorCode AdaptiveKDTree::reset_tree()
463 : : {
464 : 29 : return delete_tree_sets();
465 : : }
466 : :
467 : : } // namespace moab
468 : :
469 : : #endif
|
