aabb.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2017, Association Scientifique pour la Geologie et ses
 * Applications (ASGA). All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *     * Neither the name of ASGA nor the
 *       names of its contributors may be used to endorse or promote products
 *       derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ASGA BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 *     http://www.ring-team.org
 *
 *     RING Project
 *     Ecole Nationale Superieure de Geologie - GeoRessources
 *     2 Rue du Doyen Marcel Roubault - TSA 70605
 *     54518 VANDOEUVRE-LES-NANCY
 *     FRANCE
 */

#pragma once

#include <ringmesh/basic/box.h>
#include <ringmesh/basic/common.h>

namespace RINGMesh
{
    FORWARD_DECLARATION_DIMENSION_CLASS( MeshBase );
    FORWARD_DECLARATION_DIMENSION_CLASS( LineMesh );
    FORWARD_DECLARATION_DIMENSION_CLASS( SurfaceMeshBase );
    FORWARD_DECLARATION_DIMENSION_CLASS( VolumeMesh );
} // namespace RINGMesh

namespace RINGMesh
{
    template < index_t DIMENSION >
    class AABBTree
    {
        ringmesh_disable_copy_and_move( AABBTree );
        ringmesh_template_assert_2d_or_3d( DIMENSION );

    public:
        static const index_t ROOT_INDEX = 1;

        index_t nb_bboxes() const
        {
            return static_cast< index_t >( mapping_morton_.size() );
        }

        template < typename EvalDistance >
        std::tuple< index_t, vecn< DIMENSION >, double > closest_element_box(
            const vecn< DIMENSION >& query, const EvalDistance& action ) const
        {
            index_t nearest_box = NO_ID;
            vecn< DIMENSION > nearest_point;
            double distance;
            std::tie( nearest_box, nearest_point, distance ) =
                get_nearest_element_box_hint( query );
            closest_element_box_recursive< EvalDistance >( query, nearest_box,
                nearest_point, distance, ROOT_INDEX, 0, nb_bboxes(), action );
            ringmesh_assert( nearest_box != NO_ID );
            return std::make_tuple( nearest_box, nearest_point, distance );
        }
        /*
         * @brief Computes the intersections between a given
         * box and the element boxes.
         * @param[in] box the box to test
         * @param[in] action The functor to run when an element box intersects
         * \p box
         * @tparam EvalIntersection this functor should have an operator()
         * defined like this:
         * void operator()( index_t cur_box ) ;
         * where cur_box is the element box index
         * (e.g. in the case of AABBTree2D, this index is a polygon index)
         */
        template < class EvalIntersection >
        void compute_bbox_element_bbox_intersections(
            const Box< DIMENSION >& box, EvalIntersection& action ) const
        {
            bbox_intersect_recursive< EvalIntersection >(
                box, ROOT_INDEX, 0, nb_bboxes(), action );
        }
        /*
         * @brief Computes the self intersections of the element boxes.
         * @param[in] action The functor to run when two boxes intersect
         * @tparam EvalIntersection this functor should have an operator()
         * defined like this:
         * void operator()( index_t box1, index_t box2 ) ;
         * where box1 and box2 are the element box indices
         * (e.g. in the case of AABBTree2D, this index is a polygon index)
         */
        template < class EvalIntersection >
        void compute_self_element_bbox_intersections(
            EvalIntersection& action ) const
        {
            self_intersect_recursive< EvalIntersection >( ROOT_INDEX, 0,
                nb_bboxes(), ROOT_INDEX, 0, nb_bboxes(), action );
        }

    protected:
        AABBTree() = default;
        virtual ~AABBTree() = default;

        void initialize_tree( const std::vector< Box< DIMENSION > >& bboxes );

        bool is_leaf( index_t box_begin, index_t box_end ) const
        {
            return box_begin + 1 == box_end;
        }
        void get_recursive_iterators( index_t node_index,
            index_t box_begin,
            index_t box_end,
            index_t& middle_box,
            index_t& child_left,
            index_t& child_right ) const
        {
            middle_box = box_begin + ( box_end - box_begin ) / 2;
            child_left = 2 * node_index;
            child_right = 2 * node_index + 1;
        }

        const Box< DIMENSION >& node( index_t i ) const
        {
            ringmesh_assert( i < tree_.size() );
            return tree_[i];
        }

        Box< DIMENSION >& node( index_t i )
        {
            ringmesh_assert( i < tree_.size() );
            return tree_[i];
        }

    private:
        index_t max_node_index(
            index_t node_index, index_t box_begin, index_t box_end );
        void initialize_tree_recursive(
            const std::vector< Box< DIMENSION > >& bboxes,
            index_t node_index,
            index_t element_begin,
            index_t element_end );

        template < typename ACTION >
        void closest_element_box_recursive( const vecn< DIMENSION >& query,
            index_t& nearest_box,
            vecn< DIMENSION >& nearest_point,
            double& distance,
            index_t node_index,
            index_t element_begin,
            index_t element_end,
            const ACTION& action ) const;

        template < class ACTION >
        void bbox_intersect_recursive( const Box< DIMENSION >& box,
            index_t node_index,
            index_t element_begin,
            index_t element_end,
            ACTION& action ) const;

        template < class ACTION >
        void self_intersect_recursive( index_t node_index1,
            index_t element_begin1,
            index_t element_end1,
            index_t node_index2,
            index_t element_begin2,
            index_t element_end2,
            ACTION& action ) const;

        std::tuple< index_t, vecn< DIMENSION >, double >
            get_nearest_element_box_hint(
                const vecn< DIMENSION >& query ) const;
        virtual vecn< DIMENSION > get_point_hint_from_box(
            const Box< DIMENSION >& box, index_t element_id ) const = 0;

    protected:
        std::vector< Box< DIMENSION > > tree_{};
        std::vector< index_t > mapping_morton_{};
    };

    template < index_t DIMENSION >
    class BoxAABBTree : public AABBTree< DIMENSION >
    {
    public:
        explicit BoxAABBTree( const std::vector< Box< DIMENSION > >& boxes );

    private:
        vecn< DIMENSION > get_point_hint_from_box(
            const Box< DIMENSION >& box, index_t element_id ) const override;
    };

    ALIAS_2D_AND_3D( BoxAABBTree );

    template < index_t DIMENSION >
    class LineAABBTree : public AABBTree< DIMENSION >
    {
    public:
        explicit LineAABBTree( const LineMesh< DIMENSION >& mesh );

        std::tuple< index_t, vecn< DIMENSION >, double > closest_edge(
            const vecn< DIMENSION >& query ) const;

    private:
        vecn< DIMENSION > get_point_hint_from_box(
            const Box< DIMENSION >& box, index_t element_id ) const override;
        class DistanceToEdge
        {
        public:
            explicit DistanceToEdge( const LineMesh< DIMENSION >& mesh )
                : mesh_( mesh )
            {
            }

            std::tuple< double, vecn< DIMENSION > > operator()(
                const vecn< DIMENSION >& query, index_t cur_box ) const;

        private:
            const LineMesh< DIMENSION >& mesh_;
        };

    private:
        const LineMesh< DIMENSION >& mesh_;
    };

    ALIAS_2D_AND_3D( LineAABBTree );

    template < index_t DIMENSION >
    class SurfaceAABBTree : public AABBTree< DIMENSION >
    {
    public:
        explicit SurfaceAABBTree( const SurfaceMeshBase< DIMENSION >& mesh );

        std::tuple< index_t, vecn< DIMENSION >, double > closest_triangle(
            const vecn< DIMENSION >& query ) const;

    private:
        vecn< DIMENSION > get_point_hint_from_box(
            const Box< DIMENSION >& box, index_t element_id ) const override;
        class DistanceToTriangle
        {
        public:
            explicit DistanceToTriangle(
                const SurfaceMeshBase< DIMENSION >& mesh )
                : mesh_( mesh )
            {
            }

            std::tuple< double, vecn< DIMENSION > > operator()(
                const vecn< DIMENSION >& query, index_t cur_box ) const;

        private:
            const SurfaceMeshBase< DIMENSION >& mesh_;
        };

    private:
        const SurfaceMeshBase< DIMENSION >& mesh_;
    };

    ALIAS_2D_AND_3D( SurfaceAABBTree );

    template < index_t DIMENSION >
    class VolumeAABBTree : public AABBTree< DIMENSION >
    {
        ringmesh_template_assert_3d( DIMENSION );

    public:
        explicit VolumeAABBTree( const VolumeMesh< DIMENSION >& mesh );

        index_t containing_cell( const vecn< DIMENSION >& query ) const;

    private:
        vecn< DIMENSION > get_point_hint_from_box(
            const Box< DIMENSION >& box, index_t element_id ) const override;
        index_t containing_cell_recursive( const vecn< DIMENSION >& query,
            index_t node_index,
            index_t box_begin,
            index_t box_end ) const;

    private:
        const VolumeMesh< DIMENSION >& mesh_;
    };

    using VolumeAABBTree3D = VolumeAABBTree< 3 >;

    template < index_t DIMENSION >
    double inner_point_box_distance(
        const vecn< DIMENSION >& p, const Box< DIMENSION >& B );

    template < index_t DIMENSION >
    double point_box_signed_distance(
        const vecn< DIMENSION >& p, const Box< DIMENSION >& B );

    template < index_t DIMENSION >
    template < typename ACTION >
    void AABBTree< DIMENSION >::closest_element_box_recursive(
        const vecn< DIMENSION >& query,
        index_t& nearest_box,
        vecn< DIMENSION >& nearest_point,
        double& distance,
        index_t node_index,
        index_t box_begin,
        index_t box_end,
        const ACTION& action ) const
    {
        ringmesh_assert( node_index < tree_.size() );
        ringmesh_assert( box_begin != box_end );

        // If node is a leaf: compute point-element distance
        // and replace current if nearer
        if( is_leaf( box_begin, box_end ) )
        {
            index_t cur_box = mapping_morton_[box_begin];
            vecn< DIMENSION > cur_nearest_point;
            double cur_distance;
            std::tie( cur_distance, cur_nearest_point ) =
                action( query, cur_box );
            if( cur_distance < distance )
            {
                nearest_box = cur_box;
                nearest_point = cur_nearest_point;
                distance = cur_distance;
            }
            return;
        }
        index_t box_middle, child_left, child_right;
        get_recursive_iterators( node_index, box_begin, box_end, box_middle,
            child_left, child_right );

        double distance_left =
            point_box_signed_distance( query, node( child_left ) );
        double distance_right =
            point_box_signed_distance( query, node( child_right ) );

        // Traverse the "nearest" child first, so that it has more chances
        // to prune the traversal of the other child.
        if( distance_left < distance_right )
        {
            if( distance_left < distance )
            {
                closest_element_box_recursive< ACTION >( query, nearest_box,
                    nearest_point, distance, child_left, box_begin, box_middle,
                    action );
            }
            if( distance_right < distance )
            {
                closest_element_box_recursive< ACTION >( query, nearest_box,
                    nearest_point, distance, child_right, box_middle, box_end,
                    action );
            }
        }
        else
        {
            if( distance_right < distance )
            {
                closest_element_box_recursive< ACTION >( query, nearest_box,
                    nearest_point, distance, child_right, box_middle, box_end,
                    action );
            }
            if( distance_left < distance )
            {
                closest_element_box_recursive< ACTION >( query, nearest_box,
                    nearest_point, distance, child_left, box_begin, box_middle,
                    action );
            }
        }
    }

    template < index_t DIMENSION >
    template < typename ACTION >
    void AABBTree< DIMENSION >::bbox_intersect_recursive(
        const Box< DIMENSION >& box,
        index_t node_index,
        index_t element_begin,
        index_t element_end,
        ACTION& action ) const
    {
        ringmesh_assert( node_index < tree_.size() );
        ringmesh_assert( element_begin != element_end );

        // Prune sub-tree that does not have intersection
        if( !box.bboxes_overlap( node( node_index ) ) )
        {
            return;
        }

        // Leaf case
        if( is_leaf( element_begin, element_end ) )
        {
            // @todo Check if the box is not intersecting itself
            index_t cur_box = mapping_morton_[element_begin];
            action( cur_box );
            return;
        }

        index_t box_middle, child_left, child_right;
        get_recursive_iterators( node_index, element_begin, element_end,
            box_middle, child_left, child_right );

        bbox_intersect_recursive< ACTION >(
            box, child_left, element_begin, box_middle, action );
        bbox_intersect_recursive< ACTION >(
            box, child_right, box_middle, element_end, action );
    }

    template < index_t DIMENSION >
    template < class ACTION >
    void AABBTree< DIMENSION >::self_intersect_recursive( index_t node_index1,
        index_t element_begin1,
        index_t element_end1,
        index_t node_index2,
        index_t element_begin2,
        index_t element_end2,
        ACTION& action ) const
    {
        ringmesh_assert( element_end1 != element_begin1 );
        ringmesh_assert( element_end2 != element_begin2 );

        // Since we are intersecting the AABBTree with *itself*,
        // we can prune half of the cases by skipping the test
        // whenever node2's polygon index interval is greated than
        // node1's polygon index interval.
        if( element_end2 <= element_begin1 )
        {
            return;
        }

        // The acceleration is here:
        if( !node( node_index1 ).bboxes_overlap( node( node_index2 ) ) )
        {
            return;
        }

        // Simple case: leaf - leaf intersection.
        if( is_leaf( element_begin1, element_end1 )
            && is_leaf( element_begin2, element_end2 ) )
        {
            if( node_index1 == node_index2 )
            {
                return;
            }
            action( mapping_morton_[element_begin1],
                mapping_morton_[element_begin2] );
            return;
        }

        // If node2 has more polygons than node1, then
        //   intersect node2's two children with node1
        // else
        //   intersect node1's two children with node2
        if( element_end2 - element_begin2 > element_end1 - element_begin1 )
        {
            index_t middle_box2, child_left2, child_right2;
            get_recursive_iterators( node_index2, element_begin2, element_end2,
                middle_box2, child_left2, child_right2 );
            self_intersect_recursive< ACTION >( node_index1, element_begin1,
                element_end1, child_left2, element_begin2, middle_box2,
                action );
            self_intersect_recursive< ACTION >( node_index1, element_begin1,
                element_end1, child_right2, middle_box2, element_end2, action );
        }
        else
        {
            index_t middle_box1, child_left1, child_right1;
            get_recursive_iterators( node_index1, element_begin1, element_end1,
                middle_box1, child_left1, child_right1 );
            self_intersect_recursive< ACTION >( child_left1, element_begin1,
                middle_box1, node_index2, element_begin2, element_end2,
                action );
            self_intersect_recursive< ACTION >( child_right1, middle_box1,
                element_end1, node_index2, element_begin2, element_end2,
                action );
        }
    }
} // namespace RINGMesh