1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
/*
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

/**\file OrientedBoxTreeTool.hpp
 *\author Jason Kraftcheck ([email protected])
 *\date 2006-07-18
 */

#ifndef MOAB_ORIENTED_BOX_TREE_TOOL_HPP
#define MOAB_ORIENTED_BOX_TREE_TOOL_HPP

#include "moab/Forward.hpp"
#include "moab/GeomUtil.hpp"
#include "moab/OrientedBox.hpp"
#include <iosfwd>
#include <list>
#include <vector>

namespace moab
{

class Range;
class OrientedBox;
class StatData;
class CartVect;

/** \class OrientedBoxTreeTool
 * \brief Class for constructing and querying Hierarchical Oriented Bounding Box trees
 */
class OrientedBoxTreeTool
{
  public:
    /**\brief This provides the search range for ray intersections, measured
       relative to the origin of the ray.

       first: nonnegative limit for search
       second: negative limit for search

       These are const double* so that the window is always defined by
       pointing to other quantities, but it is not posible to change those
       quantities via the window.
     **/
    typedef std::pair< const double*, const double* > IntersectSearchWindow;

    /**\brief Misc. knobs controlling tree subdivision
     *
     * Available settings for controlling when and how nodes in the tree
     * are split.  The constructor will initialize to the default
     * settings.  All settings except best_split_ratio control when
     * a node is subdivided.  best_split_ratio influences the choice
     * of how the node is subdivided.
     *
     * A calculated ratio is used in the determination of when and how
     * to split a node.  The ratio is calculated as:
     * - \f$max(\frac{|n_L - n_R|}{n_L+n_R}, f*\frac{n_I}{n_L+n_R})\f$
     * - \f$n_L\f$ : num entities to be placed in left child
     * - \f$n_R\f$ : num entities to be placed in right child
     * - \f$f\f$ : Settings::intersect_ratio_factor
     * - \f$n_I\f$: num entities intersecting split plane
     *
     * ALL of the following conditions must be met for a node to be further
     * subdivied:
     *  - Depth must be less than max_depth
     *  - Node must contain more than max_leaf_entities entities.
     *  - The 'ratio' must be less than worst_split_ratio
     *
     * The node will be subdivided using a plane normal to one of the
     * box axis and containing the box center.  The planes are tested
     * beginning with the one orthogonal to the longest box axis and
     * finishing with the one orthogonal to the shortest box axis.  The
     * search will stop at the first plane for which the 'ratio' is
     * at least Settings::best_split_ratio .  Giving Settings::best_split_ratio
     * a non-zero value gives preference to a split orthogonal to larger
     * box dimensions.
     */
    struct Settings
    {
      public:
        Settings();             //!< set defaults
        int max_leaf_entities;  //!< Average number of entities per leaf
        int max_depth;          //!< Maximum tree depth - 0->no limit
        //! Must be in [best_split_ratio,1.0]
        //! A tree node will not be split if the ratio of children
        //! in the child nodes is greater than this value.
        double worst_split_ratio;
        //! Must be in [0.0,worst_split_ratio]
        //! The search for an optimal split plane for splitting a node
        //! will stop if at least this ratio is achieved for the number of
        //! entities on each side of the split plane.
        double best_split_ratio;
        //! Flags used to create entity sets representing tree nodes
        unsigned int set_options;
        //! Check if settings are valid.
        bool valid() const;
    };

    OrientedBoxTreeTool( Interface* i, const char* tag_name = 0, bool destroy_created_trees = false );

    ~OrientedBoxTreeTool();

    /**\brief Build oriented bounding box tree
     *
     * Build an oriented bounding box tree.
     *\param entities A list of either vertices or 2-D elements (not both)
     *                for which to build a tree.
     *\param set_handle_out A handle for the entity set representing the
     *                root of the tree.
     */
    ErrorCode build( const Range& entities, EntityHandle& set_handle_out, const Settings* settings = 0 );

    /**\brief Build a tree of sets, where each set contains triangles.
     *
     * Build a tree of sets.  Each set must contain at least one triangle
     * to define its geometry.  Each passed set will become a leaf of
     * the OBB tree.  Settings controlling tree depth are ignored by
     * this method.  The tree will be as deep as it needs to be for each
     * input set to be a leaf.
     *
     * To build a tree representing the surfaces of a geometric volume,
     * 1) Build an OBB tree for each surface using the 'build' method
     * 2) Add each surface to the contents of the resulting OBB tree root set
     * 3) Build a tree from all the surface OBB tree root sets using this
     *    method to get a combined tree for the volume.
     */
    ErrorCode join_trees( const Range& tree_roots, EntityHandle& root_set_out, const Settings* settings = 0 );

    /**\brief Traversal statistics structure
     *
     * Structure to accumulate statistics on traversal performance. Passed optionally
     * to query functions, this structure contains the count of nodes visited at each
     * level in a tree, and the count of traversals that ended at each level.
     * One TrvStats structure can be used with multiple OBB trees or multiple queries,
     * or used on only a single tree or a single query.
     *
     * Note that these traversal statistics are not related to the stats() query below,
     * which calculates static information about a tree.  These statistics relate
     * to a tree's dynamic behavior on particular operations.
     */
    class TrvStats
    {
      public:
        //! return counts of nodes visited, indexed by tree depth.
        //! the counts include both leaves and interior nodes
        const std::vector< unsigned >& nodes_visited() const
        {
            return nodes_visited_count;
        }
        //! return counts of tree leaves visited, indexed by tree depth
        const std::vector< unsigned >& leaves_visited() const
        {
            return leaves_visited_count;
        }
        //! return counts of traversals ended, indexed by tree depth
        const std::vector< unsigned >& traversals_ended() const
        {
            return traversals_ended_count;
        }
        //! return total number of ray-triangle intersection tests performed
        //! in calls made with this TrvStats
        unsigned int ray_tri_tests() const
        {
            return ray_tri_tests_count;
        }
        //! reset all counters on this structure
        void reset();
        //! print the contents of this structure to given stream
        void print( std::ostream& str ) const;

        TrvStats() : ray_tri_tests_count( 0 ) {}

      private:
        std::vector< unsigned > nodes_visited_count;
        std::vector< unsigned > leaves_visited_count;
        std::vector< unsigned > traversals_ended_count;
        unsigned int ray_tri_tests_count;

        void increment( unsigned depth );
        void increment_leaf( unsigned depth );
        void end_traversal( unsigned depth );

        friend class OrientedBoxTreeTool;
    };

    /**\brief Default/Base class to provide a context for registering intersections
     *
     * To enable different logic for how individual intersections are
     * accumulated, depending on the usage of ray_intersect_sets().
     *
     * The API to this context has 3 parts:
     * * getDesiredOrient() during initialization of ray_intersect_sets to
          determine whether this context filters by context
     * * update_orient() updates the context to know the orientation of the
         current surface wrt to its volume during a traversal visit()
     * * register_intersection() offers an intersection to the context so that
         it can decide whether to accumulate it or ignore it
     *
     * This implementation also provides a default NOP version that accumulates
     * all intersections without logic.
     *
     * A reference implementation can be found in GeomQueryTool::GQT_IntRegCtxt.
     *
     */

    class IntRegCtxt
    {

      protected:
        std::vector< double > intersections;
        std::vector< EntityHandle > sets;
        std::vector< EntityHandle > facets;

      public:
        /* provide a default behavior that will simply add the intersection data to the relevent
           lists with no logic or discrimination */
        virtual ErrorCode register_intersection( EntityHandle set,
                                                 EntityHandle tri,
                                                 double dist,
                                                 IntersectSearchWindow& /* search_win */,
                                                 GeomUtil::intersection_type /* int_type */ )
        {
            intersections.push_back( dist );
            sets.push_back( set );
            facets.push_back( tri );

            return MB_SUCCESS;
        };

        /* determine the orientation of the topological set with respect to any
           topological set in the context */
        virtual ErrorCode update_orient( EntityHandle /* set */, int* /* surfTriOrient */ )
        {
            return MB_SUCCESS;
        };

        /* determine whether or not a preferred orientation is established  */
        virtual const int* getDesiredOrient()
        {
            return NULL;
        };

        std::vector< double > get_intersections()
        {
            return intersections;
        };
        std::vector< EntityHandle > get_facets()
        {
            return facets;
        };
        std::vector< EntityHandle > get_sets()
        {
            return sets;
        };
    };

    /**\brief Intersect a ray with the triangles contained within the tree
     *
     * Intersect a ray with the triangles contained in the tree and return
     * the distance at which the intersection occured.
     *\param distances_out The output list of intersection points on the ray.
     *\param facets_out    Handles of intersected triangles corresponding to distances_out
     *\param root_set      The MBENTITYSET representing the root of the tree.
     *\param tolerance     The tolerance to use in intersection checks.
     *\param ray_point     The base point of the ray.
     *\param unit_ray_dir  The ray direction vector (must be unit length)
     *\param ray_length    Optional ray length (intersect segment instead of ray.)
     */
    ErrorCode ray_intersect_triangles( std::vector< double >& distances_out,
                                       std::vector< EntityHandle >& facets_out,
                                       EntityHandle root_set,
                                       double tolerance,
                                       const double ray_point[3],
                                       const double unit_ray_dir[3],
                                       const double* ray_length = 0,
                                       TrvStats* accum          = 0 );

    /**\brief Intersect ray with tree
     *
     * Return the tree nodes (as MBENTITYSET handles) for the leaf boxes
     * of the tree intersected by a ray.
     *\param boxes_out    The boxes intersected by the ray.
     *\param tolerance     The tolerance to use in intersection checks.
     *\param ray_point     The base point of the ray.
     *\param unit_ray_dir  The ray direction vector (must be unit length)
     *\param ray_length    Optional ray length (intersect segment instead of ray.)
     */
    ErrorCode ray_intersect_boxes( Range& boxes_out,
                                   EntityHandle root_set,
                                   double tolerance,
                                   const double ray_point[3],
                                   const double unit_ray_dir[3],
                                   const double* ray_length = 0,
                                   TrvStats* accum          = 0 );

    /**\brief Intersect ray with triangles contained in passed MBENTITYSETs
     *
     * \param raytri_test_count    If non-NULL, count of ray-triangle intersect tests
     *                             will be added to the value at which this points.
     */
    ErrorCode ray_intersect_triangles( std::vector< double >& intersection_distances_out,
                                       std::vector< EntityHandle >& intersection_facets_out,
                                       const Range& leaf_boxes_containing_tris,
                                       double tolerance,
                                       const double ray_point[3],
                                       const double unit_ray_dir[3],
                                       const double* ray_length        = 0,
                                       unsigned int* raytri_test_count = 0 );

    /**\brief Intersect a ray with the triangles contained within the tree
     *
     * Intersect a ray with the triangles contained in the tree and return
     * the distance at which the intersection occured.
     *\param distances_out The output list of intersection points on the ray.
     *\param sets_out      The contained set encountered during the tree traversal
     *                     (see 'set_build').  For the most common use, this is the
     *                     set corresponding to the geometric surface containing the
     *                     intersected triangle.
     *\param facets_out    Handles of intersected triangles corresponding to distances_out
     *\param root_set      The MBENTITYSET representing the root of the tree.
     *\param tolerance     The tolerance to use in intersection checks.
     *\param ray_point     The base point of the ray.
     *\param unit_ray_dir  The ray direction vector (must be unit length)
     *\param search_win    An interval that defines the current window in which the
     *                     an intersection is being sought: (nonnegative, negative)
     *\param register_intersection A context for assessing and registering intersections
     *                     derived from IntRegCtxt
     *\param accum         Optional class for tree traversal statistics.
     */

    ErrorCode ray_intersect_sets( std::vector< double >& distances_out,
                                  std::vector< EntityHandle >& sets_out,
                                  std::vector< EntityHandle >& facets_out,
                                  EntityHandle root_set,
                                  double tolerance,
                                  const double ray_point[3],
                                  const double unit_ray_dir[3],
                                  IntersectSearchWindow& search_win,
                                  IntRegCtxt& register_intersection,
                                  TrvStats* accum = 0 );

    /*\brief Version that doesn't require a search window or an intersection registration context
     *
     *
     */
    ErrorCode ray_intersect_sets( std::vector< double >& distances_out,
                                  std::vector< EntityHandle >& sets_out,
                                  std::vector< EntityHandle >& facets_out,
                                  EntityHandle root_set,
                                  double tolerance,
                                  const double ray_point[3],
                                  const double unit_ray_dir[3],
                                  const double* ray_length = 0,
                                  TrvStats* accum          = 0 );

    ErrorCode ray_intersect_sets( EntityHandle root_set,
                                  double tolerance,
                                  const double ray_point[3],
                                  const double unit_ray_dir[3],
                                  IntersectSearchWindow& search_win,
                                  IntRegCtxt& register_intersection,
                                  TrvStats* accum = 0 );

    /**\brief Find closest surface, facet in surface, and location on facet
     *
     * Find the closest location in the tree to the specified location.
     *\param point Location to search from
     *\param point_out Closest location on closest facet
     *\param facet_out Closest 2D element to input position
     *\param set_out Set containing closest facet.  0 if tree was not
     *               constructed using 'set_build'
     */
    ErrorCode closest_to_location( const double* point,
                                   EntityHandle tree_root,
                                   double* point_out,
                                   EntityHandle& facet_out,
                                   EntityHandle* set_out = 0,
                                   TrvStats* accum       = 0 );

    /**\brief Find closest facet(s) to input position.
     *
     * Find the closest location(s) in the tree to the specified location.
     *\param point Location to search from
     *\param facets_out Closest 2D elements to input position are appended to this list
     *\param sets_out If non-null, sets owning facets are appended to this list.
     */
    ErrorCode closest_to_location( const double* point,
                                   EntityHandle tree_root,
                                   double tolerance,
                                   std::vector< EntityHandle >& facets_out,
                                   std::vector< EntityHandle >* sets_out = 0,
                                   TrvStats* accum                       = 0 );

    /**\brief Find facets intersected by a sphere
     *
     * Find facets intersected by a spherical volume.
     *\param center     Sphere center
     *\param radius     Sphere radius
     *\param tree_root  Root of OBB tree
     *\param facets_out List of triangles intersecting sphere
     *\param sets_out   If not null, sets owning facets are appended to this
     *                  list in an order corresponding to the entries in
     *                  facets_out.
     */
    ErrorCode sphere_intersect_triangles( const double* center,
                                          double radius,
                                          EntityHandle tree_root,
                                          std::vector< EntityHandle >& facets_out,
                                          std::vector< EntityHandle >* sets_out = 0,
                                          TrvStats* accum                       = 0 );

    ErrorCode get_close_tris( CartVect int_pt,
                              double tol,
                              const EntityHandle* rootSet,
                              const EntityHandle* geomVol,
                              const Tag* senseTag,
                              std::vector< EntityHandle >& close_tris,
                              std::vector< int >& close_senses );

    /**\brief Get oriented box at node in tree
     *
     * Get the oriented box for a node in an oriented bounding box tree.
     */
    ErrorCode box( EntityHandle node_set, double center[3], double axis1[3], double axis2[3], double axis3[3] );

    ErrorCode delete_tree( EntityHandle root_set );

    /**\brief Remove obb tree root from the Oriented Box Tree Tool data structure
     *
     * Remove obb tree root from the Oriented Box Tree Tool data structure (createdTrees)
     */
    ErrorCode remove_root( EntityHandle root_set );

    /**\brief Print out tree
     *
     * Print the tree to an output stream in a human-readable form.
     *\param tree_root_set  Entity set representing tree root.
     *\param list_contents  If true, list entities in each tree node,
     *                      If false, just list number of entities.
     *\param id_tag_name    If specified, must be the name of an existing
     *                      integer tag containing an ID for the entities.
     *                      Not used if list_contents is false.
     */
    void print( EntityHandle tree_root_set,
                std::ostream& stream,
                bool list_contents      = false,
                const char* id_tag_name = 0 );

    /**\brief Print tree statistics
     *
     * Print misc. stats. describing tree
     */
    ErrorCode stats( EntityHandle tree_root_set, std::ostream& stream );

    /**\brief Get tree statistics
     *
     * Get summary stats. describing tree
     * \param set Root of tree for which data is requested
     * \param total_entities Entities in tree
     * \param root_volume Total volume of root box
     * \param tot_node_volume Total volume in all nodes
     * \param tot_to_root_volume Ratio of total / root volume
     * \param tree_height Maximum height of tree, from root to leaf
     * \param node_count Number of nodes in tree
     * \param num_leaves Number of leaf nodes in tree
     */
    ErrorCode stats( EntityHandle set,
                     unsigned& entities_in_tree,
                     double& root_volume,
                     double& tot_node_volume,
                     double& tot_to_root_volume,
                     unsigned& tree_height,
                     unsigned& node_count,
                     unsigned& num_leaves );

    /** \brief Implement this and pass instance to preorder_traverse
     *
     * This interface may be implemented and an instance passed to
     * preorder_traverse to define some operation to do when traversing
     * the tree.
     */
    class Op
    {
      public:
        /**\brief Visit a node in the tree during a traversal.
         *
         * This method is called for each node in the tree visited
         * during a pre-order traversal.
         *\param node The EntityHandle for the entity set for the tree node.
         *\param depth The current depth in the tree.
         *\param descend Output: if false, traversal will skip children
         *             of the current node, or if the current node is a
         *             leaf, the 'leaf' method will not be called.
         */
        virtual ErrorCode visit( EntityHandle node, int depth, bool& descend ) = 0;

        /**\brief Process a leaf node during tree traversal */
        virtual ErrorCode leaf( EntityHandle node ) = 0;

        virtual ~Op();  // probably isn't necessary in this case, and
                        // does nothing, but lots of compilers warn if
                        // virtual function but no virtual destructor.
    };

    /**\brief Visitor pattern - do operation for each tree node
     *
     * Do a preorder traversal of the tree, calling the method
     * in the passed operation instance for each node in the tree.
     * Parent node is visited before either child (pre-order traversal).
     * If operator method passes back the 'descend' argument as false,
     * traversal will not descend to the children of the current node.
     */
    ErrorCode preorder_traverse( EntityHandle root_set, Op& operation, TrvStats* accum = 0 );

    Interface* get_moab_instance() const
    {
        return instance;
    }

    struct SetData;

    /**\brief Get oriented box at node in tree
     *
     * Get the oriented box for a node in an oriented bounding box tree.
     *
     * NOTE: This function is provided for internal MOAB use only.
     *       The definition of OrientedBox is not available as a part
     *       of the MOAB API
     */
    ErrorCode box( EntityHandle node_set, OrientedBox& box );

  private:
    ErrorCode build_tree( const Range& entities, EntityHandle& set, int depth, const Settings& settings );

    ErrorCode build_sets( std::list< SetData >& sets, EntityHandle& node_set, int depth, const Settings& settings );

    ErrorCode recursive_stats( OrientedBoxTreeTool* tool,
                               Interface* instance,
                               EntityHandle set,
                               int depth,
                               StatData& data,
                               unsigned& count_out,
                               CartVect& dimensions_out );

    Interface* instance;
    Tag tagHandle;

    bool cleanUpTrees;
    std::vector< EntityHandle > createdTrees;
};

}  // namespace moab

#endif