Branch data Line data Source code
1 : : /**\file BVHTree.hpp
2 : : * \class moab::BVHTree
3 : : * \brief Bounding Volume Hierarchy (sorta like a tree), for sorting and searching entities
4 : : * spatially
5 : : */
6 : :
7 : : #ifndef BVH_TREE_HPP
8 : : #define BVH_TREE_HPP
9 : :
10 : : #include "moab/Interface.hpp"
11 : : #include "moab/CartVect.hpp"
12 : : #include "moab/BoundBox.hpp"
13 : : #include "moab/Tree.hpp"
14 : : #include "moab/Range.hpp"
15 : : #include "moab/CN.hpp"
16 : :
17 : : #include <vector>
18 : : #include <cfloat>
19 : : #include <climits>
20 : : #include <map>
21 : : #include <set>
22 : : #include <iostream>
23 : :
24 : : #include "moab/win32_config.h"
25 : :
26 : : namespace moab
27 : : {
28 : :
29 : : class ElemEvaluator;
30 : :
31 : : class BVHTree : public Tree
32 : : {
33 : : public:
34 : : BVHTree( Interface* impl );
35 : :
36 : 1 : ~BVHTree()
37 : 1 : {
38 : 1 : reset_tree();
39 [ - + ]: 1 : }
40 : :
41 : : /** \brief Destroy the tree maintained by this object, optionally checking we have the right
42 : : * root. \param root If non-NULL, check that this is the root, return failure if not
43 : : */
44 : : virtual ErrorCode reset_tree();
45 : :
46 : : virtual ErrorCode parse_options( FileOptions& opts );
47 : :
48 : : /** \brief Get bounding box for tree below tree_node, or entire tree
49 : : * If no tree has been built yet, returns +/- DBL_MAX for all dimensions. Note for some tree
50 : : * types, boxes are not available for non-root nodes, and this function will return failure if
51 : : * non-root is passed in \param box The box for this tree \param tree_node If non-NULL, node for
52 : : * which box is requested, tree root if NULL \return Only returns error on fatal condition
53 : : */
54 : : virtual ErrorCode get_bounding_box( BoundBox& box, EntityHandle* tree_node = NULL ) const;
55 : :
56 : : /** Build the tree
57 : : * Build a tree with the entities input. If a non-NULL tree_root_set pointer is input,
58 : : * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct
59 : : * a new root set and pass its handle back in *tree_root_set. Options vary by tree type;
60 : : * see Tree.hpp for common options; options specific to AdaptiveKDTree:
61 : : * SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
62 : : * CANDIDATE_PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
63 : : * for possible values; default = 1 (SUBDIVISION_SNAP)
64 : : * \param entities Entities with which to build the tree
65 : : * \param tree_root Root set for tree (see function description)
66 : : * \param opts Options for tree (see function description)
67 : : * \return Error is returned only on build failure
68 : : */
69 : : virtual ErrorCode build_tree( const Range& entities, EntityHandle* tree_root_set = NULL,
70 : : FileOptions* options = NULL );
71 : :
72 : : /** \brief Get leaf containing input position.
73 : : *
74 : : * Does not take into account global bounding box of tree.
75 : : * - Therefore there is always one leaf containing the point.
76 : : * - If caller wants to account for global bounding box, then
77 : : * caller can test against that box and not call this method
78 : : * at all if the point is outside the box, as there is no leaf
79 : : * containing the point in that case.
80 : : * \param point Point to be located in tree
81 : : * \param leaf_out Leaf containing point
82 : : * \param iter_tol Tolerance for convergence of point search
83 : : * \param inside_tol Tolerance for inside element calculation
84 : : * \param multiple_leaves Some tree types can have multiple leaves containing a point;
85 : : * if non-NULL, this parameter is returned true if multiple leaves contain
86 : : * the input point
87 : : * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
88 : : * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
89 : : */
90 : : virtual ErrorCode point_search( const double* point, EntityHandle& leaf_out, const double iter_tol = 1.0e-10,
91 : : const double inside_tol = 1.0e-6, bool* multiple_leaves = NULL,
92 : : EntityHandle* start_node = NULL, CartVect* params = NULL );
93 : :
94 : : /** \brief Find all leaves within a given distance from point
95 : : * If dists_out input non-NULL, also returns distances from each leaf; if
96 : : * point i is inside leaf, 0 is given as dists_out[i].
97 : : * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities
98 : : * in tree nodes and return containing entities in leaves_out. In those cases, if params_out
99 : : * is also non-NULL, will return parameters in those elements in that vector.
100 : : * \param from_point Point to be located in tree
101 : : * \param distance Distance within which to query
102 : : * \param result_list Leaves within distance or containing point
103 : : * \param iter_tol Tolerance for convergence of point search
104 : : * \param inside_tol Tolerance for inside element calculation
105 : : * \param result_dists If non-NULL, will contain distsances to leaves
106 : : * \param result_params If non-NULL, will contain parameters of the point in the ents in
107 : : * leaves_out \param tree_root Start from this tree node (non-NULL) instead of tree root (NULL)
108 : : */
109 : : virtual ErrorCode distance_search( const double from_point[3], const double distance,
110 : : std::vector< EntityHandle >& result_list, const double iter_tol = 1.0e-10,
111 : : const double inside_tol = 1.0e-6, std::vector< double >* result_dists = NULL,
112 : : std::vector< CartVect >* result_params = NULL, EntityHandle* tree_root = NULL );
113 : :
114 : : //! print various things about this tree
115 : : virtual ErrorCode print();
116 : :
117 : : private:
118 : : // don't allow copy constructor, too complicated
119 : : BVHTree( const BVHTree& s );
120 : :
121 : 84678 : class HandleData
122 : : {
123 : : public:
124 : 729 : HandleData( EntityHandle h, const BoundBox& bx, const double dp ) : myHandle( h ), myBox( bx ), myDim( dp ) {}
125 : : HandleData() : myHandle( 0 ), myDim( -1 ) {}
126 : : EntityHandle myHandle;
127 : : BoundBox myBox;
128 : : double myDim;
129 : : };
130 : : typedef std::vector< HandleData > HandleDataVec;
131 : :
132 : 7144 : class SplitData
133 : : {
134 : : public:
135 : 188 : SplitData()
136 [ + - ][ + - ]: 188 : : dim( UINT_MAX ), nl( UINT_MAX ), nr( UINT_MAX ), split( DBL_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX )
137 : : {
138 : 188 : }
139 : 3384 : SplitData( const SplitData& f )
140 : : : dim( f.dim ), nl( f.nl ), nr( f.nr ), split( f.split ), Lmax( f.Lmax ), Rmin( f.Rmin ),
141 : 3384 : boundingBox( f.boundingBox ), leftBox( f.leftBox ), rightBox( f.rightBox )
142 : : {
143 : 3384 : }
144 : : unsigned int dim, nl, nr;
145 : : double split;
146 : : double Lmax, Rmin;
147 : : BoundBox boundingBox, leftBox, rightBox;
148 : 188 : SplitData& operator=( const SplitData& f )
149 : : {
150 : 188 : dim = f.dim;
151 : 188 : nl = f.nl;
152 : 188 : nr = f.nr;
153 : 188 : split = f.split;
154 : 188 : Lmax = f.Lmax;
155 : 188 : Rmin = f.Rmin;
156 : 188 : boundingBox = f.boundingBox;
157 : 188 : leftBox = f.leftBox;
158 : 188 : rightBox = f.rightBox;
159 : 188 : return *this;
160 : : }
161 : : }; // SplitData
162 : :
163 : : class Split_comparator : public std::binary_function< SplitData, SplitData, bool >
164 : : {
165 : 3008 : inline double objective( const SplitData& a ) const
166 : : {
167 : 3008 : return a.Lmax * a.nl - a.Rmin * a.nr;
168 : : }
169 : :
170 : : public:
171 : 1504 : bool operator()( const SplitData& a, const SplitData& b ) const
172 : : {
173 : 1504 : return objective( a ) < objective( b );
174 : : }
175 : : }; // Split_comparator
176 : :
177 : : class HandleData_comparator : public std::binary_function< HandleData, HandleData, bool >
178 : : {
179 : : public:
180 : 46484 : bool operator()( const HandleData& a, const HandleData& b )
181 : : {
182 [ + + ][ + + ]: 46484 : return a.myDim < b.myDim || ( a.myDim == b.myDim && a.myHandle < b.myHandle );
[ + + ]
183 : : }
184 : : }; // Iterator_comparator
185 : :
186 : 6016 : class Bucket
187 : : {
188 : : public:
189 : 1504 : Bucket() : mySize( 0 ) {}
190 : 2256 : Bucket( const Bucket& f ) : mySize( f.mySize ), boundingBox( f.boundingBox ) {}
191 : : Bucket( const unsigned int sz ) : mySize( sz ) {}
192 : : static unsigned int bucket_index( int num_splits, const BoundBox& box, const BoundBox& interval,
193 : : const unsigned int dim );
194 : : unsigned int mySize;
195 : : BoundBox boundingBox;
196 : : Bucket& operator=( const Bucket& f )
197 : : {
198 : : boundingBox = f.boundingBox;
199 : : mySize = f.mySize;
200 : : return *this;
201 : : }
202 : : }; // Bucket
203 : :
204 : 3714 : class Node
205 : : {
206 : : public:
207 : : HandleDataVec entities;
208 : : unsigned int dim, child;
209 : : double Lmax, Rmin;
210 : : BoundBox box;
211 [ + - ]: 377 : Node() : dim( UINT_MAX ), child( UINT_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX ) {}
212 : : Node& operator=( const Node& f )
213 : : {
214 : : dim = f.dim;
215 : : child = f.child;
216 : : Lmax = f.Lmax;
217 : : Rmin = f.Rmin;
218 : : entities = f.entities;
219 : : return *this;
220 : : }
221 : : }; // Node
222 : :
223 : 1508 : class TreeNode
224 : : {
225 : : public:
226 : : unsigned int dim, child;
227 : : double Lmax, Rmin;
228 : : BoundBox box;
229 : 377 : TreeNode( int dm, int chld, double lmx, double rmn, BoundBox& bx )
230 : 377 : : dim( dm ), child( chld ), Lmax( lmx ), Rmin( rmn ), box( bx )
231 : : {
232 : 377 : }
233 : : TreeNode& operator=( const TreeNode& f )
234 : : {
235 : : dim = f.dim;
236 : : child = f.child;
237 : : Lmax = f.Lmax;
238 : : Rmin = f.Rmin;
239 : : return *this;
240 : : }
241 : : }; // TreeNode
242 : :
243 : : void establish_buckets( HandleDataVec::const_iterator begin, HandleDataVec::const_iterator end,
244 : : const BoundBox& interval, std::vector< std::vector< Bucket > >& buckets ) const;
245 : :
246 : : unsigned int set_interval( BoundBox& interval, std::vector< Bucket >::const_iterator begin,
247 : : std::vector< Bucket >::const_iterator end ) const;
248 : :
249 : : void initialize_splits( std::vector< std::vector< SplitData > >& splits,
250 : : const std::vector< std::vector< Bucket > >& buckets, const SplitData& data ) const;
251 : :
252 : : void order_elements( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
253 : :
254 : : void median_order( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
255 : :
256 : : void choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const;
257 : :
258 : : void find_split( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
259 : :
260 : : ErrorCode find_point( const std::vector< double >& point, const unsigned int& index, const double iter_tol,
261 : : const double inside_tol, std::pair< EntityHandle, CartVect >& result );
262 : :
263 : : EntityHandle bruteforce_find( const double* point, const double iter_tol = 1.0e-10,
264 : : const double inside_tol = 1.0e-6 );
265 : :
266 : : int local_build_tree( std::vector< Node >& tree_nodes, HandleDataVec::iterator begin, HandleDataVec::iterator end,
267 : : const int index, const BoundBox& box, const int depth = 0 );
268 : :
269 : : // builds up vector of HandleData, which caches elements' bounding boxes
270 : : ErrorCode construct_element_vec( std::vector< HandleData >& handle_data_vec, const Range& elements,
271 : : BoundBox& bounding_box );
272 : :
273 : : // convert the std::vector<Node> to myTree and a bunch of entity sets
274 : : ErrorCode convert_tree( std::vector< Node >& tree_nodes );
275 : :
276 : : // print tree nodes
277 : : ErrorCode print_nodes( std::vector< Node >& nodes );
278 : :
279 : : Range entityHandles;
280 : : std::vector< TreeNode > myTree;
281 : : int splitsPerDir;
282 : : EntityHandle startSetHandle;
283 : : static MOAB_EXPORT const char* treeName;
284 : : }; // class Bvh_tree
285 : :
286 : 44052 : inline unsigned int BVHTree::Bucket::bucket_index( int num_splits, const BoundBox& box, const BoundBox& interval,
287 : : const unsigned int dim )
288 : : {
289 : : // see FastMemoryEfficientCellLocationinUnstructuredGridsForVisualization.pdf
290 : : // around page 9
291 : :
292 : : // Paper arithmetic is over-optimized.. this is safer.
293 : 44052 : const double min = interval.bMin[dim];
294 : 44052 : const double length = ( interval.bMax[dim] - min ) / ( num_splits + 1 );
295 : 44052 : const double center = ( ( box.bMax[dim] + box.bMin[dim] ) / 2.0 ) - min;
296 : : #ifndef NDEBUG
297 : : #ifdef BVH_SHOW_INDEX
298 : : std::cout << "[" << min << " , " << interval.max[dim] << " ]" << std::endl;
299 : : std::cout << "[" << box.bMin[dim] << " , " << box.bMax[dim] << " ]" << std::endl;
300 : : std::cout << "Length of bucket" << length << std::endl;
301 : : std::cout << "Center: " << ( box.bMax[dim] + box.bMin[dim] ) / 2.0 << std::endl;
302 : : std::cout << "Distance of center from min: " << center << std::endl;
303 : : std::cout << "ratio: " << center / length << std::endl;
304 : : std::cout << "index: " << std::ceil( center / length ) - 1 << std::endl;
305 : : #endif
306 : : #endif
307 : 44052 : unsigned int cl = std::ceil( center / length );
308 [ + - ]: 44052 : return ( cl > 0 ? cl - 1 : 0 );
309 : : }
310 : :
311 [ + - ][ + - ]: 2 : inline BVHTree::BVHTree( Interface* impl ) : Tree( impl ), splitsPerDir( 3 ), startSetHandle( 0 )
312 : : {
313 [ + - ]: 1 : boxTagName = treeName;
314 : 1 : }
315 : :
316 : 3384 : inline unsigned int BVHTree::set_interval( BoundBox& interval, std::vector< Bucket >::const_iterator begin,
317 : : std::vector< Bucket >::const_iterator end ) const
318 : : {
319 : 3384 : bool first = true;
320 : 3384 : unsigned int count = 0;
321 [ + - ][ + - ]: 10152 : for( std::vector< Bucket >::const_iterator b = begin; b != end; ++b )
[ + + ]
322 : : {
323 [ + - ]: 6768 : const BoundBox& box = b->boundingBox;
324 [ + - ]: 6768 : count += b->mySize;
325 [ + - ][ + + ]: 6768 : if( b->mySize != 0 )
326 : : {
327 [ + + ]: 4227 : if( first )
328 : : {
329 [ + - ]: 2811 : interval = box;
330 : 2811 : first = false;
331 : : }
332 : : else
333 [ + - ]: 4227 : interval.update( box );
334 : : }
335 : : }
336 : 3384 : return count;
337 : : }
338 : :
339 : 96 : inline void BVHTree::order_elements( HandleDataVec::iterator& begin, HandleDataVec::iterator& end,
340 : : SplitData& data ) const
341 : : {
342 [ + - ][ + - ]: 5160 : for( HandleDataVec::iterator i = begin; i != end; ++i )
[ + + ]
343 : : {
344 [ + - ][ + - ]: 5064 : const int index = Bucket::bucket_index( splitsPerDir, i->myBox, data.boundingBox, data.dim );
345 [ + - ]: 5064 : i->myDim = ( index <= data.split ) ? 0 : 1;
346 : : }
347 [ + - ]: 96 : std::sort( begin, end, HandleData_comparator() );
348 : 96 : }
349 : :
350 : 188 : inline void BVHTree::choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const
351 : : {
352 [ + - ]: 188 : std::vector< SplitData > best_splits;
353 [ + - ][ + - ]: 752 : for( std::vector< std::vector< SplitData > >::const_iterator i = splits.begin(); i != splits.end(); ++i )
[ + + ]
354 : : {
355 : : std::vector< SplitData >::const_iterator j =
356 [ + - ][ + - ]: 564 : std::min_element( ( *i ).begin(), ( *i ).end(), Split_comparator() );
[ + - ]
357 [ + - ][ + - ]: 564 : best_splits.push_back( *j );
358 : : }
359 [ + - ][ + - ]: 188 : data = *std::min_element( best_splits.begin(), best_splits.end(), Split_comparator() );
[ + - ]
360 : 188 : }
361 : :
362 : 1 : inline ErrorCode BVHTree::construct_element_vec( std::vector< HandleData >& handle_data_vec, const Range& elements,
363 : : BoundBox& bounding_box )
364 : : {
365 [ + - ]: 1 : std::vector< double > coordinate( 3 * CN::MAX_NODES_PER_ELEMENT );
366 : : int num_conn;
367 : 1 : ErrorCode rval = MB_SUCCESS;
368 [ + - ]: 2 : std::vector< EntityHandle > storage;
369 [ + - ][ + - ]: 730 : for( Range::iterator i = elements.begin(); i != elements.end(); ++i )
[ + - ][ + - ]
[ + + ]
370 : : {
371 : : // TODO: not generic enough. Why dim != 3
372 : : // Commence un-necessary deep copying.
373 : : const EntityHandle* connect;
374 [ + - ][ + - ]: 729 : rval = mbImpl->get_connectivity( *i, connect, num_conn, false, &storage );
375 [ - + ]: 729 : if( MB_SUCCESS != rval ) return rval;
376 [ + - ][ + - ]: 729 : rval = mbImpl->get_coords( connect, num_conn, &coordinate[0] );
377 [ - + ]: 729 : if( MB_SUCCESS != rval ) return rval;
378 [ + - ]: 729 : BoundBox box;
379 [ + + ]: 6561 : for( int j = 0; j < num_conn; j++ )
380 [ + - ][ + - ]: 5832 : box.update( &coordinate[3 * j] );
381 [ + - ][ + - ]: 729 : if( i == elements.begin() )
[ + + ]
382 [ + - ]: 1 : bounding_box = box;
383 : : else
384 [ + - ]: 728 : bounding_box.update( box );
385 [ + - ][ + - ]: 729 : handle_data_vec.push_back( HandleData( *i, box, 0.0 ) );
[ + - ]
386 : 729 : }
387 : :
388 : 2 : return rval;
389 : : }
390 : :
391 : 1 : inline ErrorCode BVHTree::reset_tree()
392 : : {
393 : 1 : return delete_tree_sets();
394 : : }
395 : :
396 : : } // namespace moab
397 : :
398 : : #endif // BVH_TREE_HPP
|
