mesh_AABB.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2022 Inria
 *  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 the ALICE Project-Team 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 THE COPYRIGHT HOLDER OR CONTRIBUTORS 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.
 *
 *  Contact: Bruno Levy
 *
 *     https://www.inria.fr/fr/bruno-levy
 *
 *     Inria,
 *     Domaine de Voluceau,
 *     78150 Le Chesnay - Rocquencourt
 *     FRANCE
 *
 */
 
#ifndef GEOGRAM_MESH_MESH_AABB
#define GEOGRAM_MESH_MESH_AABB
 
#include <geogram/basic/common.h>
#include <geogram/mesh/mesh.h>
#include <geogram/basic/geometry.h>
#include <geogram/basic/stopwatch.h>
 
namespace GEO {
 
    enum AABBReorderMode {
        AABB_INPLACE, AABB_INDIRECT, AABB_NOREORDER
    };
 
    template <class BOX> class AABB {
 
    public:
        void enlarge_boxes(double amount) {
            for(auto& B: bboxes_) {
                B.enlarge(amount);
            }
        }
 
    protected:
 
        void initialize(
            index_t nb, std::function<void(BOX&, index_t)> get_bbox
        ) {
            nb_ = nb;
            bboxes_.resize(max_node_index(1, 0, nb) + 1);
            // +1 because size == max_index + 1 !!!
            init_bboxes_recursive(1, 0, nb_, get_bbox);
        }
 
        void bbox_intersect_recursive(
            std::function<void(index_t)> action, const BOX& box,
            index_t node, index_t b, index_t e
        ) const {
            geo_debug_assert(e != b);
 
            // Prune sub-tree that does not have intersection
            if(!bboxes_overlap(box, bboxes_[node])) {
                return;
            }
 
            // Leaf case
            if(e == b+1) {
                action(element_in_leaf(b));
                return;
            }
 
            // Recursion
            index_t m = b + (e - b) / 2;
            index_t node_l = 2 * node;
            index_t node_r = 2 * node + 1;
 
            bbox_intersect_recursive(action, box, node_l, b, m);
            bbox_intersect_recursive(action, box, node_r, m, e);
        }
 
        void self_intersect_recursive(
            std::function<void(index_t,index_t)> action,
            index_t node1, index_t b1, index_t e1,
            index_t node2, index_t b2, index_t e2
        ) const {
            geo_debug_assert(e1 != b1);
            geo_debug_assert(e2 != b2);
 
            // Since we are intersecting the AABBTree with *itself*,
            // we can prune half of the cases by skipping the test
            // whenever node2's facet index interval is greater than
            // node1's facet index interval.
            if(e2 <= b1) {
                return;
            }
 
            // The acceleration is here:
            if(
                (node1 != node2) &&
                !bboxes_overlap(bboxes_[node1], bboxes_[node2])
            ) {
                return;
            }
 
            // Simple case: leaf - leaf intersection.
            if(b1 + 1 == e1 && b2 + 1 == e2) {
                if(b1 != b2) {
                    action(element_in_leaf(b1), element_in_leaf(b2));
                }
                return;
            }
 
            // If node2 has more elements than node1, then
            //   intersect node2's two children with node1
            // else
            //   intersect node1's two children with node2
            if(e2 - b2 > e1 - b1) {
                index_t m2 = b2 + (e2 - b2) / 2;
                index_t node2_l = 2 * node2;
                index_t node2_r = 2 * node2 + 1;
                self_intersect_recursive(action, node1, b1, e1, node2_l, b2, m2);
                self_intersect_recursive(action, node1, b1, e1, node2_r, m2, e2);
            } else {
                index_t m1 = b1 + (e1 - b1) / 2;
                index_t node1_l = 2 * node1;
                index_t node1_r = 2 * node1 + 1;
                self_intersect_recursive(action, node1_l, b1, m1, node2, b2, e2);
                self_intersect_recursive(action, node1_r, m1, e1, node2, b2, e2);
            }
        }
 
        void other_intersect_recursive(
            std::function<void(index_t,index_t)> action,
            index_t node1, index_t b1, index_t e1,
            const AABB<BOX>* other,
            index_t node2, index_t b2, index_t e2
        ) const {
            geo_debug_assert(e1 != b1);
            geo_debug_assert(e2 != b2);
 
            // The acceleration is here:
            if(!bboxes_overlap(bboxes_[node1], other->bboxes_[node2])) {
                return;
            }
 
            // Simple case: leaf - leaf intersection.
            if(b1 + 1 == e1 && b2 + 1 == e2) {
                action(element_in_leaf(b1), element_in_leaf(b2));
                return;
            }
 
            // If node2 has more elements than node1, then
            //   intersect node2's two children with node1
            // else
            //   intersect node1's two children with node2
            if(e2 - b2 > e1 - b1) {
                index_t m2 = b2 + (e2 - b2) / 2;
                index_t node2_l = 2 * node2;
                index_t node2_r = 2 * node2 + 1;
                other_intersect_recursive(
                    action, node1, b1, e1, other, node2_l, b2, m2
                );
                other_intersect_recursive(
                    action, node1, b1, e1, other, node2_r, m2, e2
                );
            } else {
                index_t m1 = b1 + (e1 - b1) / 2;
                index_t node1_l = 2 * node1;
                index_t node1_r = 2 * node1 + 1;
                other_intersect_recursive(
                    action, node1_l, b1, m1, other, node2, b2, e2
                );
                other_intersect_recursive(
                    action, node1_r, m1, e1, other, node2, b2, e2
                );
            }
        }
 
 
        static index_t max_node_index(index_t node_index, index_t b, index_t e) {
            geo_debug_assert(e > b);
            if(b + 1 == e) {
                return node_index;
            }
            index_t m = b + (e - b) / 2;
            index_t childl = 2 * node_index;
            index_t childr = 2 * node_index + 1;
            return std::max(
                max_node_index(childl, b, m),
                max_node_index(childr, m, e)
            );
        }
 
        void init_bboxes_recursive(
            index_t node_index,
            index_t b, index_t e,
            std::function<void(BOX&, index_t)> get_bbox
        ) {
            geo_debug_assert(node_index < bboxes_.size());
            geo_debug_assert(b != e);
            if(b + 1 == e) {
                get_bbox(bboxes_[node_index], element_in_leaf(b));
                return;
            }
            index_t m = b + (e - b) / 2;
            index_t childl = 2 * node_index;
            index_t childr = 2 * node_index + 1;
            geo_debug_assert(childl < bboxes_.size());
            geo_debug_assert(childr < bboxes_.size());
            init_bboxes_recursive(childl, b, m, get_bbox);
            init_bboxes_recursive(childr, m, e, get_bbox);
            geo_debug_assert(childl < bboxes_.size());
            geo_debug_assert(childr < bboxes_.size());
            bbox_union(bboxes_[node_index], bboxes_[childl], bboxes_[childr]);
        }
 
        bool indirect() const {
            return (reorder_.size() != 0);
        }
 
        index_t element_in_leaf(index_t i) const {
            return indirect() ? reorder_[i] : i;
        }
 
    protected:
        index_t nb_;
        vector<BOX> bboxes_;
        vector<index_t> reorder_; 
    };
 
    typedef AABB<Box2d> AABB2d;
    typedef AABB<Box3d> AABB3d;
 
    /**************************************************************/
 
    class GEOGRAM_API MeshAABB2d : public AABB2d {
    public:
        MeshAABB2d() : mesh_(nullptr) {
        }
 
        const Mesh* mesh() const {
            return mesh_;
        }
 
    protected:
        Mesh* mesh_;
    };
 
    /**************************************************************/
 
    class GEOGRAM_API MeshAABB3d : public AABB3d {
    public:
        MeshAABB3d() : mesh_(nullptr) {
        }
 
        const Mesh* mesh() const {
            return mesh_;
        }
 
    protected:
        Mesh* mesh_;
    };
 
    /**************************************************************/
 
    class GEOGRAM_API MeshFacetsAABB : public MeshAABB3d {
    public:
 
        struct Intersection {
            Intersection() :
                t(Numeric::max_float64()),
                f(NO_INDEX),
                i(NO_INDEX), j(NO_INDEX), k(NO_INDEX)
                {
                }
            vec3 p;        
            double t;      
            index_t f;     
            vec3 N;        
            index_t i,j,k; 
            double u,v;    
        };
 
 
        MeshFacetsAABB() {
        }
 
        MeshFacetsAABB(Mesh& M, AABBReorderMode reorder_mode) {
            initialize(M, reorder_mode);
        }
 
        MeshFacetsAABB(const Mesh& M) {
            initialize(const_cast<Mesh&>(M), AABB_INDIRECT);
        }
 
        void initialize(Mesh& M, AABBReorderMode reorder_mode = AABB_INDIRECT);
 
#ifndef GOMGEN
        [[deprecated("use MeshFacetsAABB(Mesh&,AABBReorderMode) instead")]]
        MeshFacetsAABB(Mesh& M, bool reorder) {
            initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
        }
 
        [[deprecated("use initialize(Mesh&,AABBReorderMode) instead")]]
        void initialize(Mesh& M, bool reorder) {
            initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
        }
#endif
        void compute_facet_bbox_intersections(
            std::function<void(index_t, index_t)> action,
            bool concurrent = false
        ) const {
            if(
                Process::maximum_concurrent_threads() <= 1 ||
                mesh_->facets.nb() <= 1024
            ) {
                self_intersect_recursive(
                    action,
                    1, 0, mesh_->facets.nb(),
                    1, 0, mesh_->facets.nb()
                );
            } else {
                self_bbox_intersections_parallel(action, concurrent);
            }
 
        }
 
        void compute_bbox_facet_bbox_intersections(
            const Box& box_in, std::function<void(index_t)> action
        ) const {
            bbox_intersect_recursive(
                action, box_in, 1, 0, mesh_->facets.nb()
            );
        }
 
        index_t nearest_facet(
            const vec3& p, vec3& nearest_point, double& sq_dist
        ) const {
            if(mesh_->facets.nb() == 0) {
                sq_dist = Numeric::max_float64();
                return NO_INDEX;
            }
            index_t nearest_facet;
            get_nearest_facet_hint(p, nearest_facet, nearest_point, sq_dist);
            nearest_facet_recursive(
                p,
                nearest_facet, nearest_point, sq_dist,
                1, 0, mesh_->facets.nb()
            );
            return nearest_facet;
        }
 
        index_t nearest_facet(const vec3& p) const {
            vec3 nearest_point;
            double sq_dist;
            return nearest_facet(p, nearest_point, sq_dist);
        }
 
 
        index_t nearest_facet_filtered(
            const vec3& p, vec3& nearest_point, double& sq_dist,
            std::function<bool(index_t)> filter
        ) const {
            if(mesh_->facets.nb() == 0) {
                return NO_INDEX;
            }
            index_t nearest_facet = NO_INDEX;
            sq_dist = Numeric::max_float64();
            nearest_facet_recursive_filtered(
                p,
                nearest_facet, nearest_point, sq_dist,
                1, 0, mesh_->facets.nb(),
                filter
            );
            return nearest_facet;
        }
 
        void nearest_facet_with_hint(
            const vec3& p,
            index_t& nearest_facet, vec3& nearest_point, double& sq_dist
        ) const {
            if(mesh_->facets.nb() == 0) {
                nearest_facet = NO_INDEX;
                sq_dist = Numeric::max_float64();
                nearest_point=vec3(0,0,0);
                return;
            }
            if(nearest_facet == NO_FACET) {
                get_nearest_facet_hint(
                    p, nearest_facet, nearest_point, sq_dist
                );
            }
            nearest_facet_recursive(
                p,
                nearest_facet, nearest_point, sq_dist,
                1, 0, mesh_->facets.nb()
            );
        }
 
        double squared_distance(const vec3& p) const {
            vec3 nearest_point;
            double result;
            nearest_facet(p, nearest_point, result);
            return result;
        }
 
        bool ray_intersection(
            const Ray& R,
            double tmax = Numeric::max_float64(),
            index_t ignore_f = NO_INDEX
        ) const;
 
 
        bool ray_nearest_intersection(const Ray& R, Intersection& I) const;
 
        bool segment_intersection(const vec3& q1, const vec3& q2) const {
            return ray_intersection(Ray(q1, q2-q1), 1.0);
        }
 
        bool segment_nearest_intersection(
            const vec3& q1, const vec3& q2, double& t, index_t& f
        ) const {
            Ray R(q1, q2-q1);
            Intersection I;
            I.t = 1.0;
            bool result = ray_nearest_intersection(R,I);
            t = I.t;
            f = I.f;
            return result;
        }
 
        void ray_all_intersections(
            const Ray& R,
            std::function<void(const Intersection&)> action
        ) const;
 
 
        void line_all_intersections(
            const vec3& O, const vec3& D,
            std::function<void(const Intersection&)> action
        ) const;
 
 
        bool contains(const vec3& p) const;
 
    protected:
 
        void get_nearest_facet_hint(
            const vec3& p,
            index_t& nearest_facet, vec3& nearest_point, double& sq_dist
        ) const;
 
        void nearest_facet_recursive(
            const vec3& p,
            index_t& nearest_facet, vec3& nearest_point, double& sq_dist,
            index_t n, index_t b, index_t e
        ) const;
 
 
        void nearest_facet_recursive_filtered(
            const vec3& p,
            index_t& nearest_facet, vec3& nearest_point, double& sq_dist,
            index_t n, index_t b, index_t e,
            std::function<bool(index_t)> filter
        ) const;
 
        bool ray_intersection_recursive(
            const Ray& R, const vec3& dirinv, double max_t, index_t ignore_f,
            index_t n, index_t b, index_t e
        ) const;
 
        void ray_nearest_intersection_recursive(
            const Ray& R, const vec3& dirinv, Intersection& I, index_t ignore_f,
            index_t n, index_t b, index_t e, index_t coord
        ) const;
 
 
        void ray_all_intersections_recursive(
            const Ray& R, const vec3& dirinv,
            std::function<void(const Intersection&)> action,
            index_t n, index_t b, index_t e
        ) const;
 
 
        void line_all_intersections_recursive(
            const vec3& O, const vec3& D, const vec3& dirinv,
            std::function<void(const Intersection&)> action,
            index_t n, index_t b, index_t e
        ) const;
 
 
        void self_bbox_intersections_parallel(
            std::function<void(index_t, index_t)> action, bool concurrent=false
        ) const;
    };
 
    /***********************************************************************/
 
    class GEOGRAM_API MeshCellsAABB : public MeshAABB3d {
    public:
 
        static constexpr index_t NO_TET = NO_INDEX;
 
        MeshCellsAABB() {
        }
 
        MeshCellsAABB(Mesh& M, AABBReorderMode reorder_mode) {
            initialize(M, reorder_mode);
        }
 
        MeshCellsAABB(const Mesh& M) {
            initialize(const_cast<Mesh&>(M), AABB_INDIRECT);
        }
 
        void initialize(Mesh& M, AABBReorderMode reorder_mode = AABB_INDIRECT);
 
#ifndef GOMGEN
        [[deprecated("use MeshCellsAABB(Mesh&,AABBReorderMode) instead")]]
        MeshCellsAABB(Mesh& M, bool reorder) {
            initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
        }
 
        [[deprecated("use initialize(Mesh&,AABBReorderMode) instead")]]
        void initialize(Mesh& M, bool reorder) {
            initialize(M, reorder ? AABB_INPLACE : AABB_NOREORDER);
        }
#endif
 
        index_t containing_tet(const vec3& p) const {
            geo_debug_assert(mesh_->cells.are_simplices());
            return containing_tet_recursive(
                p, 1, 0, mesh_->cells.nb()
            );
        }
 
        void compute_bbox_cell_bbox_intersections(
            const Box& box_in,
            std::function<void(index_t)> action
        ) const {
            bbox_intersect_recursive(
                action, box_in, 1, 0, mesh_->cells.nb()
            );
        }
 
        void containing_boxes(
            const vec3& p, std::function<void(index_t)> action
        ) const {
            containing_bboxes_recursive(
                action, p, 1, 0, mesh_->cells.nb()
            );
        }
 
        void compute_cell_bbox_intersections(
            std::function<void(index_t, index_t)> action
        ) const {
            self_intersect_recursive(
                action,
                1, 0, mesh_->cells.nb(),
                1, 0, mesh_->cells.nb()
            );
        }
 
        void compute_other_cell_bbox_intersections(
            MeshCellsAABB* other,
            std::function<void(index_t, index_t)> action
        ) const {
            other_intersect_recursive(
                action,
                1, 0, mesh_->cells.nb(),
                other,
                1, 0, other->mesh_->cells.nb()
            );
        }
 
 
    protected:
 
        index_t containing_tet_recursive(
            const vec3& p,
            index_t n, index_t b, index_t e
        ) const;
 
 
        void containing_bboxes_recursive(
            std::function<void(index_t)> action,
            const vec3& p,
            index_t node, index_t b, index_t e
        ) const {
            geo_debug_assert(e != b);
 
            // Prune sub-tree that does not have intersection
            if(!bboxes_[node].contains(p)) {
                return;
            }
 
            // Leaf case
            if(e == b+1) {
                action(element_in_leaf(b));
                return;
            }
 
            // Recursion
            index_t m = b + (e - b) / 2;
            index_t node_l = 2 * node;
            index_t node_r = 2 * node + 1;
 
            containing_bboxes_recursive(action, p, node_l, b, m);
            containing_bboxes_recursive(action, p, node_r, m, e);
        }
    };
 
/*******************************************************************/
 
    class GEOGRAM_API MeshFacetsAABB2d : public MeshAABB2d {
    public:
 
        static constexpr index_t NO_TRIANGLE = NO_INDEX;
 
        MeshFacetsAABB2d();
 
        MeshFacetsAABB2d(Mesh& M, bool reorder = true);
 
        void initialize(Mesh& M, bool reorder = true);
 
        index_t containing_triangle(const vec2& p) const {
            geo_debug_assert(mesh_->facets.are_simplices());
            return containing_triangle_recursive(
                p, 1, 0, mesh_->facets.nb()
            );
        }
 
        void compute_bbox_cell_bbox_intersections(
            const Box2d& box_in,
            std::function<void(index_t)> action
        ) const {
            bbox_intersect_recursive(
                action, box_in, 1, 0, mesh_->facets.nb()
            );
        }
 
        void containing_boxes(
            const vec2& p, std::function<void(index_t)> action
        ) const {
            containing_bboxes_recursive(
                action, p, 1, 0, mesh_->facets.nb()
            );
        }
 
        void compute_facet_bbox_intersections(
            std::function<void(index_t, index_t)> action
        ) const {
            self_intersect_recursive(
                action,
                1, 0, mesh_->facets.nb(),
                1, 0, mesh_->facets.nb()
            );
        }
 
        void compute_other_cell_bbox_intersections(
            MeshFacetsAABB2d* other,
            std::function<void(index_t, index_t)> action
        ) const {
            other_intersect_recursive(
                action,
                1, 0, mesh_->facets.nb(),
                other,
                1, 0, other->mesh_->facets.nb()
            );
        }
 
 
    protected:
 
        index_t containing_triangle_recursive(
            const vec2& p,
            index_t n, index_t b, index_t e
        ) const;
 
 
        void containing_bboxes_recursive(
            std::function<void(index_t)> action,
            const vec2& p,
            index_t node, index_t b, index_t e
        ) const {
            geo_debug_assert(e != b);
 
            // Prune sub-tree that does not have intersection
            if(!bboxes_[node].contains(p)) {
                return;
            }
 
            // Leaf case
            if(e == b+1) {
                action(element_in_leaf(b));
                return;
            }
 
            // Recursion
            index_t m = b + (e - b) / 2;
            index_t node_l = 2 * node;
            index_t node_r = 2 * node + 1;
 
            containing_bboxes_recursive(action, p, node_l, b, m);
            containing_bboxes_recursive(action, p, node_r, m, e);
        }
    };
 
 
}
 
#endif