OrientedBoxTreeTool.hpp

Go to the documentation of this file.

/*
 * 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 (kraftche@cae.wisc.edu)
 *\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