Cppcheck report - MOAB: src/moab/BVHTree.hpp

/**\file BVHTree.hpp
 * \class moab::BVHTree
 * \brief Bounding Volume Hierarchy (sorta like a tree), for sorting and searching entities
 * spatially
 */

#ifndef BVH_TREE_HPP
#define BVH_TREE_HPP

#include "moab/Interface.hpp"
#include "moab/CartVect.hpp"
#include "moab/BoundBox.hpp"
#include "moab/Tree.hpp"
#include "moab/Range.hpp"
#include "moab/CN.hpp"

#include <vector>
#include <cfloat>
#include <climits>
#include <map>
#include <set>
#include <iostream>

#include "moab/win32_config.h"

namespace moab
{

class ElemEvaluator;

class BVHTree : public Tree
{
  public:
    BVHTree( Interface* impl );
<--- Class 'BVHTree' has a constructor with 1 argument that is not explicit. [+]
Class 'BVHTree' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
<--- Class 'BVHTree' has a constructor with 1 argument that is not explicit. [+]
Class 'BVHTree' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.

    ~BVHTree()
    {
        reset_tree();
    }

    /** \brief Destroy the tree maintained by this object, optionally checking we have the right
     * root. \param root If non-NULL, check that this is the root, return failure if not
     */
    virtual ErrorCode reset_tree();

    virtual ErrorCode parse_options( FileOptions& opts );

    /** \brief Get bounding box for tree below tree_node, or entire tree
     * If no tree has been built yet, returns +/- DBL_MAX for all dimensions.  Note for some tree
     * types, boxes are not available for non-root nodes, and this function will return failure if
     * non-root is passed in \param box The box for this tree \param tree_node If non-NULL, node for
     * which box is requested, tree root if NULL \return Only returns error on fatal condition
     */
    virtual ErrorCode get_bounding_box( BoundBox& box, EntityHandle* tree_node = NULL ) const;

    /** Build the tree
     * Build a tree with the entities input.  If a non-NULL tree_root_set pointer is input,
     * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct
     * a new root set and pass its handle back in *tree_root_set.  Options vary by tree type;
     * see Tree.hpp for common options; options specific to AdaptiveKDTree:
     * SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
     * CANDIDATE_PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
     *          for possible values; default = 1 (SUBDIVISION_SNAP)
     * \param entities Entities with which to build the tree
     * \param tree_root Root set for tree (see function description)
     * \param opts Options for tree (see function description)
     * \return Error is returned only on build failure
     */
    virtual ErrorCode build_tree( const Range& entities,
                                  EntityHandle* tree_root_set = NULL,
                                  FileOptions* options        = NULL );

    /** \brief Get leaf containing input position.
     *
     * Does not take into account global bounding box of tree.
     * - Therefore there is always one leaf containing the point.
     * - If caller wants to account for global bounding box, then
     * caller can test against that box and not call this method
     * at all if the point is outside the box, as there is no leaf
     * containing the point in that case.
     * \param point Point to be located in tree
     * \param leaf_out Leaf containing point
     * \param iter_tol Tolerance for convergence of point search
     * \param inside_tol Tolerance for inside element calculation
     * \param multiple_leaves Some tree types can have multiple leaves containing a point;
     *          if non-NULL, this parameter is returned true if multiple leaves contain
     *          the input point
     * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
     * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
     */
    virtual ErrorCode point_search( const double* point,
                                    EntityHandle& leaf_out,
                                    const double iter_tol    = 1.0e-10,
                                    const double inside_tol  = 1.0e-6,
                                    bool* multiple_leaves    = NULL,
                                    EntityHandle* start_node = NULL,
                                    CartVect* params         = NULL );

    /** \brief Find all leaves within a given distance from point
     * If dists_out input non-NULL, also returns distances from each leaf; if
     * point i is inside leaf, 0 is given as dists_out[i].
     * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities
     * in tree nodes and return containing entities in leaves_out.  In those cases, if params_out
     * is also non-NULL, will return parameters in those elements in that vector.
     * \param from_point Point to be located in tree
     * \param distance Distance within which to query
     * \param result_list Leaves within distance or containing point
     * \param iter_tol Tolerance for convergence of point search
     * \param inside_tol Tolerance for inside element calculation
     * \param result_dists If non-NULL, will contain distsances to leaves
     * \param result_params If non-NULL, will contain parameters of the point in the ents in
     * leaves_out \param tree_root Start from this tree node (non-NULL) instead of tree root (NULL)
     */
    virtual ErrorCode distance_search( const double from_point[3],
                                       const double distance,
                                       std::vector< EntityHandle >& result_list,
                                       const double iter_tol                  = 1.0e-10,
                                       const double inside_tol                = 1.0e-6,
                                       std::vector< double >* result_dists    = NULL,
                                       std::vector< CartVect >* result_params = NULL,
                                       EntityHandle* tree_root                = NULL );

    //! print various things about this tree
    virtual ErrorCode print();

  private:
    // don't allow copy constructor, too complicated
    BVHTree( const BVHTree& s );

    class HandleData
    {
      public:
        HandleData( EntityHandle h, const BoundBox& bx, const double dp ) : myHandle( h ), myBox( bx ), myDim( dp ) {}
        HandleData() : myHandle( 0 ), myDim( -1 ) {}
        EntityHandle myHandle;
        BoundBox myBox;
        double myDim;
    };
    typedef std::vector< HandleData > HandleDataVec;

    class SplitData
    {
      public:
        SplitData()
            : dim( UINT_MAX ), nl( UINT_MAX ), nr( UINT_MAX ), split( DBL_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX )
        {
        }
        SplitData( const SplitData& f )
            : dim( f.dim ), nl( f.nl ), nr( f.nr ), split( f.split ), Lmax( f.Lmax ), Rmin( f.Rmin ),
              boundingBox( f.boundingBox ), leftBox( f.leftBox ), rightBox( f.rightBox )
        {
        }
        unsigned int dim, nl, nr;
        double split;
        double Lmax, Rmin;
        BoundBox boundingBox, leftBox, rightBox;
        SplitData& operator=( const SplitData& f )
        {
            dim         = f.dim;
            nl          = f.nl;
            nr          = f.nr;
            split       = f.split;
            Lmax        = f.Lmax;
            Rmin        = f.Rmin;
            boundingBox = f.boundingBox;
            leftBox     = f.leftBox;
            rightBox    = f.rightBox;
            return *this;
        }
    };  // SplitData

    class Split_comparator
    {
        // deprecation of binary_function
        typedef SplitData first_argument_type;
        typedef SplitData second_argument_type;
        typedef bool result_type;

        inline double objective( const SplitData& a ) const
        {
            return a.Lmax * a.nl - a.Rmin * a.nr;
        }

      public:
        bool operator()( const SplitData& a, const SplitData& b ) const
        {
            return objective( a ) < objective( b );
        }
    };  // Split_comparator

    class HandleData_comparator
    {
        // deprecation of binary_function
        typedef HandleData first_argument_type;
        typedef HandleData second_argument_type;
        typedef bool result_type;

      public:
        bool operator()( const HandleData& a, const HandleData& b )
        {
            return a.myDim < b.myDim || ( a.myDim == b.myDim && a.myHandle < b.myHandle );
        }
    };  // Iterator_comparator

    class Bucket
    {
      public:
        Bucket() : mySize( 0 ) {}
        Bucket( const Bucket& f ) : mySize( f.mySize ), boundingBox( f.boundingBox ) {}
        Bucket( const unsigned int sz ) : mySize( sz ) {}
<--- Class 'Bucket' has a constructor with 1 argument that is not explicit. [+]
Class 'Bucket' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
<--- Class 'Bucket' has a constructor with 1 argument that is not explicit. [+]
Class 'Bucket' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
        static unsigned int bucket_index( int num_splits,
                                          const BoundBox& box,
                                          const BoundBox& interval,
                                          const unsigned int dim );
        unsigned int mySize;
        BoundBox boundingBox;
        Bucket& operator=( const Bucket& f )
        {
            boundingBox = f.boundingBox;
            mySize      = f.mySize;
            return *this;
        }
    };  // Bucket

    class Node
    {
      public:
        HandleDataVec entities;
        unsigned int dim, child;
        double Lmax, Rmin;
        BoundBox box;
        Node() : dim( UINT_MAX ), child( UINT_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX ) {}
        Node( const Node& f ) : entities( f.entities ), dim( f.dim ), child( f.child ), Lmax( f.Lmax ), Rmin( f.Rmin )
        {
        }
        Node& operator=( const Node& f )
        {
            dim      = f.dim;
            child    = f.child;
            Lmax     = f.Lmax;
            Rmin     = f.Rmin;
            entities = f.entities;
            return *this;
        }
    };  // Node

    class TreeNode
    {
      public:
        unsigned int dim, child;
        double Lmax, Rmin;
        BoundBox box;
        TreeNode( int dm, int chld, double lmx, double rmn, BoundBox& bx )
            : dim( dm ), child( chld ), Lmax( lmx ), Rmin( rmn ), box( bx )
        {
        }
        TreeNode( const TreeNode& f )
        {
            dim   = f.dim;
            child = f.child;
            Lmax  = f.Lmax;
            Rmin  = f.Rmin;
        }
        TreeNode& operator=( const TreeNode& f )
        {
            dim   = f.dim;
            child = f.child;
            Lmax  = f.Lmax;
            Rmin  = f.Rmin;
            return *this;
        }
    };  // TreeNode

    void establish_buckets( HandleDataVec::const_iterator begin,
                            HandleDataVec::const_iterator end,
                            const BoundBox& interval,
                            std::vector< std::vector< Bucket > >& buckets ) const;

    unsigned int set_interval( BoundBox& interval,
                               std::vector< Bucket >::const_iterator begin,
                               std::vector< Bucket >::const_iterator end ) const;

    void initialize_splits( std::vector< std::vector< SplitData > >& splits,
                            const std::vector< std::vector< Bucket > >& buckets,
                            const SplitData& data ) const;

    void order_elements( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;

    void median_order( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;

    void choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const;

    void find_split( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;

    ErrorCode find_point( const std::vector< double >& point,
                          const unsigned int& index,
                          const double iter_tol,
                          const double inside_tol,
                          std::pair< EntityHandle, CartVect >& result );

    EntityHandle bruteforce_find( const double* point,
                                  const double iter_tol   = 1.0e-10,
                                  const double inside_tol = 1.0e-6 );

    int local_build_tree( std::vector< Node >& tree_nodes,
                          HandleDataVec::iterator begin,
                          HandleDataVec::iterator end,
                          const int index,
                          const BoundBox& box,
                          const int depth = 0 );

    // builds up vector of HandleData, which caches elements' bounding boxes
    ErrorCode construct_element_vec( std::vector< HandleData >& handle_data_vec,
                                     const Range& elements,
                                     BoundBox& bounding_box );

    // convert the std::vector<Node> to myTree and a bunch of entity sets
    ErrorCode convert_tree( std::vector< Node >& tree_nodes );

    // print tree nodes
    ErrorCode print_nodes( std::vector< Node >& nodes );

    Range entityHandles;
    std::vector< TreeNode > myTree;
    int splitsPerDir;
    EntityHandle startSetHandle;
    static MOAB_EXPORT const char* treeName;
};  // class Bvh_tree

inline unsigned int BVHTree::Bucket::bucket_index( int num_splits,
                                                   const BoundBox& box,
                                                   const BoundBox& interval,
                                                   const unsigned int dim )
{
    // see FastMemoryEfficientCellLocationinUnstructuredGridsForVisualization.pdf
    // around page 9

    // Paper arithmetic is over-optimized.. this is safer.
    const double min    = interval.bMin[dim];
    const double length = ( interval.bMax[dim] - min ) / ( num_splits + 1 );
    const double center = ( ( box.bMax[dim] + box.bMin[dim] ) / 2.0 ) - min;
#ifndef NDEBUG
#ifdef BVH_SHOW_INDEX
    std::cout << "[" << min << " , " << interval.max[dim] << " ]" << std::endl;
    std::cout << "[" << box.bMin[dim] << " , " << box.bMax[dim] << " ]" << std::endl;
    std::cout << "Length of bucket" << length << std::endl;
    std::cout << "Center: " << ( box.bMax[dim] + box.bMin[dim] ) / 2.0 << std::endl;
    std::cout << "Distance of center from min:  " << center << std::endl;
    std::cout << "ratio: " << center / length << std::endl;
    std::cout << "index: " << std::ceil( center / length ) - 1 << std::endl;
#endif
#endif
    unsigned int cl = std::ceil( center / length );
    return ( cl > 0 ? cl - 1 : 0 );
}

inline BVHTree::BVHTree( Interface* impl ) : Tree( impl ), splitsPerDir( 3 ), startSetHandle( 0 )
{
    boxTagName = treeName;
}

inline unsigned int BVHTree::set_interval( BoundBox& interval,
                                           std::vector< Bucket >::const_iterator begin,
                                           std::vector< Bucket >::const_iterator end ) const
{
    bool first         = true;
    unsigned int count = 0;
    for( std::vector< Bucket >::const_iterator b = begin; b != end; ++b )
    {
        const BoundBox& box = b->boundingBox;
        count += b->mySize;
        if( b->mySize != 0 )
        {
            if( first )
            {
                interval = box;
                first    = false;
            }
            else
                interval.update( box );
        }
    }
    return count;
}

inline void BVHTree::order_elements( HandleDataVec::iterator& begin,
                                     HandleDataVec::iterator& end,
                                     SplitData& data ) const
{
    for( HandleDataVec::iterator i = begin; i != end; ++i )
    {
        const int index = Bucket::bucket_index( splitsPerDir, i->myBox, data.boundingBox, data.dim );
        i->myDim        = ( index <= data.split ) ? 0 : 1;
    }
    std::sort( begin, end, HandleData_comparator() );
}

inline void BVHTree::choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const
{
    std::vector< SplitData > best_splits;
    for( std::vector< std::vector< SplitData > >::const_iterator i = splits.begin(); i != splits.end(); ++i )
    {
        std::vector< SplitData >::const_iterator j =
            std::min_element( ( *i ).begin(), ( *i ).end(), Split_comparator() );
        best_splits.push_back( *j );
    }
    data = *std::min_element( best_splits.begin(), best_splits.end(), Split_comparator() );
}

inline ErrorCode BVHTree::construct_element_vec( std::vector< HandleData >& handle_data_vec,
                                                 const Range& elements,
                                                 BoundBox& bounding_box )
{
    std::vector< double > coordinate( 3 * CN::MAX_NODES_PER_ELEMENT );
    int num_conn;
    ErrorCode rval = MB_SUCCESS;
    std::vector< EntityHandle > storage;
    for( Range::iterator i = elements.begin(); i != elements.end(); ++i )
    {
        // TODO: not generic enough. Why dim != 3
        // Commence un-necessary deep copying.
        const EntityHandle* connect;
        rval = mbImpl->get_connectivity( *i, connect, num_conn, false, &storage );
        if( MB_SUCCESS != rval ) return rval;
        rval = mbImpl->get_coords( connect, num_conn, &coordinate[0] );
        if( MB_SUCCESS != rval ) return rval;
        BoundBox box;
        for( int j = 0; j < num_conn; j++ )
            box.update( &coordinate[3 * j] );
        if( i == elements.begin() )
            bounding_box = box;
        else
            bounding_box.update( box );
        handle_data_vec.push_back( HandleData( *i, box, 0.0 ) );
    }

    return rval;
}

inline ErrorCode BVHTree::reset_tree()
{
    return delete_tree_sets();
}

}  // namespace moab

#endif  // BVH_TREE_HPP