mesh_surface_intersection.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_SURFACE_INTERSECTION
#define GEOGRAM_MESH_MESH_SURFACE_INTERSECTION
 
#include <geogram/basic/common.h>
#include <geogram/mesh/mesh.h>
#include <geogram/mesh/mesh_io.h>
#include <geogram/numerics/exact_geometry.h>
#include <geogram/numerics/predicates.h>
#include <geogram/basic/process.h>
#include <geogram/basic/attributes.h>
#include <geogram/basic/debug_stream.h>
#include <functional>
#include <tuple>
 
namespace GEO {
 
    struct IsectInfo;
 
    /********************************************************************/
 
    class GEOGRAM_API MeshSurfaceIntersection {
    public:
 
    typedef exact::vec3h ExactPoint;
 
    MeshSurfaceIntersection(Mesh& M);
    ~MeshSurfaceIntersection();
 
    void intersect();
 
    void remove_external_shell();
 
    void remove_internal_shells();
 
    void remove_fins();
 
    void classify(const std::string& expr);
 
    index_t compute_component_inclusion_bits(index_t component, index_t v);
 
    index_t compute_component_inclusion_bits_exact(index_t component, index_t v);
 
    index_t compute_component_inclusion_bits_exact_exact(
        index_t component, index_t v
    );
 
    void simplify_coplanar_facets(double angle_tolerance = 0.0);
 
    void set_verbose(bool x) {
        verbose_ = x;
        if(!verbose_ && fine_verbose_) {
            fine_verbose_ = false;
        }
    }
 
    void set_fine_verbose(bool x) {
        fine_verbose_ = x;
        if(fine_verbose_ && !verbose_) {
            verbose_ = true;
        }
    }
 
    void set_monster_threshold(index_t nb) {
        monster_threshold_ = nb;
    }
 
    void set_dry_run(bool x) {
        dry_run_ = x;
    }
 
    void set_delaunay(bool x) {
        delaunay_ = x;
    }
 
    void set_detect_intersecting_neighbors(bool x) {
        detect_intersecting_neighbors_ = x;
    }
 
    void set_radial_sort(bool x) {
        use_radial_sort_ = x;
    }
 
    void set_build_skeleton(Mesh* skeleton, bool trim_fins=false) {
        skeleton_ = skeleton;
        skeleton_trim_fins_ = trim_fins;
    }
 
    void set_interpolate_attributes(bool x) {
        interpolate_attributes_ = x;
    }
 
    protected:
    void intersect_prologue();
 
    void intersect_get_intersections(vector<IsectInfo>& intersections);
 
    void intersect_remesh_intersections(vector<IsectInfo>& intersections);
 
    void intersect_epilogue(const vector<IsectInfo>& intersections);
 
 
    void lock() {
        Process::acquire_spinlock(lock_);
    }
 
    void unlock() {
        Process::release_spinlock(lock_);
    }
 
    ExactPoint exact_vertex(index_t v) const;
 
    bool is_original_vertex(index_t v) const {
        return (vertex_to_exact_point_[v] == nullptr);
    }
 
    index_t find_or_create_exact_vertex(const ExactPoint& p);
 
    Mesh& target_mesh() {
        return mesh_;
    }
 
    const Mesh& target_mesh() const {
        return mesh_;
    }
 
    const Mesh& readonly_mesh() const {
        return mesh_copy_;
    }
 
    class RadialSort;
 
    void build_Weiler_model();
 
    void mark_external_shell(vector<index_t>& on_external_shell);
 
    index_t get_initial_facet(index_t f) const {
        return original_facet_id_[f];
    }
 
    bool initial_facet_is_flipped(index_t f) const {
        return !f_is_flipped_.is_bound() || f_is_flipped_[f];
    }
 
    std::tuple<vec3, vec3, vec3> get_initial_facet_vertices(index_t f) const {
        // All facets are duplicated before radial sort. If f is one of
        // the duplicated facets then its orientation is flipped as
        // compared to initial facet.
        index_t orig_f = original_facet_id_[f];
        vec3 p1 = mesh_copy_.facets.point(orig_f,0);
        vec3 p2 = mesh_copy_.facets.point(orig_f,1);
        vec3 p3 = mesh_copy_.facets.point(orig_f,2);
        if(initial_facet_is_flipped(f)) {
            std::swap(p1,p3);
        }
        return std::make_tuple(p1,p2,p3);
    }
 
    protected:
 
    class GEOGRAM_API RadialSort {
    public:
        RadialSort(const MeshSurfaceIntersection& I) :
            I_(I),
            mesh_(I_.target_mesh()),
            h_ref_(NO_INDEX),
            degenerate_(false)
        {
        }
 
        void init(index_t h_ref);
 
        bool operator()(index_t h1, index_t h2) const;
 
        bool degenerate() const {
            return degenerate_;
        }
 
    protected:
        Sign h_orient(index_t h1, index_t h2) const;
 
        Sign h_refNorient(index_t h2) const;
 
        exact::vec3 normal(index_t h) const;
 
        void report_problem(const char* message) const;
 
    private:
        const MeshSurfaceIntersection& I_;
        const Mesh& mesh_;
        index_t h_ref_;     // reference halfedge
        exact::vec3 N_ref_; // normal to reference triangle (exact)
        mutable bool degenerate_;
    };
 
    protected:
    Process::spinlock lock_;
    Mesh& mesh_;
    Mesh mesh_copy_;
    Attribute<const ExactPoint*> vertex_to_exact_point_;
    Attribute<index_t> original_facet_id_; // mesh_ facet to mesh_copy_ facet
    Attribute<bool> f_is_flipped_;   // mesh_ facet is flipped wrt mesh_copy_
 
#if defined(GEOGRAM_USE_EXACT_NT) && defined(GEOGRAM_EXACT_NT_IS_MPF_NT)
    // Exact points are canonicalized
    // (by Numeric::optimize_number_representation(vec3HEx)) so
    // we can use this comparator that makes the global vertex map
    // much much faster.
    typedef vec3HExLexicoCompareCanonical ExactPointCompare;
#else
    // Generic comparator for global vertex map.
    typedef vec3HgLexicoCompare<exact::scalar> ExactPointCompare;
#endif
    std::map<ExactPoint,index_t,ExactPointCompare> exact_point_to_vertex_;
 
    bool verbose_;
    bool fine_verbose_;
    bool delaunay_;
    bool detect_intersecting_neighbors_;
    bool use_radial_sort_;
 
    PCK::SOSMode SOS_bkp_;
 
    index_t monster_threshold_;
    bool dry_run_;
    friend class MeshInTriangle;
    friend class CoplanarFacets;
 
    Mesh* skeleton_;
    bool skeleton_trim_fins_;
    bool interpolate_attributes_;
 
    bool has_operand_bits_;
    /***************************************************/
 
    class Halfedges {
    public:
 
        Halfedges(MeshSurfaceIntersection& I) : mesh_(I.mesh_) {
        }
 
        ~Halfedges() {
            // TODO: destroy alpha3 attribute (kept now for debugging
        }
 
        void initialize() {
            facet_corner_alpha3_.bind(
                mesh_.facet_corners.attributes(), "alpha3"
            );
        }
 
        index_t nb() const {
            return mesh_.facet_corners.nb();
        }
 
        index_as_iterator begin() const {
            return index_as_iterator(0);
        }
 
        index_as_iterator end() const {
            return index_as_iterator(nb());
        }
 
        index_t facet(index_t h) const {
            return h/3;
        }
 
        index_t alpha2(index_t h) const {
            index_t t1 = facet(h);
            index_t t2 = mesh_.facet_corners.adjacent_facet(h);
            if(t2 == NO_INDEX) {
                return NO_INDEX;
            }
            for(index_t h2: mesh_.facets.corners(t2)) {
                if(mesh_.facet_corners.adjacent_facet(h2) == t1) {
                    return h2;
                }
            }
            geo_assert_not_reached;
        }
 
        index_t alpha3(index_t h) const {
            return facet_corner_alpha3_[h];
        }
 
        index_t facet_alpha3(index_t f) const {
            return alpha3(3*f)/3;
        }
 
        index_t vertex(index_t h, index_t dlv) const {
            index_t f  = h/3;
            index_t lv = (h+dlv)%3;
            return mesh_.facets.vertex(f,lv);
        }
 
 
        void sew2(index_t h1, index_t h2) {
            geo_debug_assert(vertex(h1,0) == vertex(h2,1));
            geo_debug_assert(vertex(h2,0) == vertex(h1,1));
            index_t t1 = h1/3;
            index_t t2 = h2/3;
            mesh_.facet_corners.set_adjacent_facet(h1,t2);
            mesh_.facet_corners.set_adjacent_facet(h2,t1);
        }
 
        void sew3(index_t h1, index_t h2) {
            geo_debug_assert(vertex(h1,0) == vertex(h2,1));
            geo_debug_assert(vertex(h2,0) == vertex(h1,1));
            facet_corner_alpha3_[h1] = h2;
            facet_corner_alpha3_[h2] = h1;
        }
 
    private:
        Mesh& mesh_;
        Attribute<index_t> facet_corner_alpha3_;
    } halfedges_;
 
    /***************************************************/
 
    class RadialBundles {
    public:
 
        RadialBundles(MeshSurfaceIntersection& I) : I_(I), mesh_(I.mesh_) {
        }
 
        void initialize();
 
        index_t nb() const {
            return bndl_start_.size() - 1;
        }
 
        index_as_iterator begin() const {
            return index_as_iterator(0);
        }
 
        index_as_iterator end() const {
            return index_as_iterator(nb());
        }
 
        index_t nb_halfedges(index_t bndl) const {
            geo_debug_assert(bndl < nb());
            return bndl_start_[bndl+1] - bndl_start_[bndl];
        }
 
        index_t halfedge(index_t bndl, index_t li) const {
            geo_debug_assert(bndl < nb());
            geo_debug_assert(li < nb_halfedges(bndl));
            return H_[bndl_start_[bndl] + li];
        }
 
        void set_halfedge(index_t bndl, index_t li, index_t h) {
            geo_debug_assert(bndl < nb());
            geo_debug_assert(li < nb_halfedges(bndl));
            H_[bndl_start_[bndl] + li] = h;
        }
 
        index_ptr_range halfedges(index_t bndl) {
            return index_ptr_range(
                H_, bndl_start_[bndl], bndl_start_[bndl+1]
            );
        }
 
        const_index_ptr_range halfedges(index_t bndl) const {
            return const_index_ptr_range(
                H_, bndl_start_[bndl], bndl_start_[bndl+1]
            );
        }
 
        index_t vertex(index_t bndl, index_t lv) const {
            geo_debug_assert(bndl_start_[bndl+1] - bndl_start_[bndl] > 0);
            index_t h = H_[bndl_start_[bndl]];
            return I_.halfedges_.vertex(h,lv);
        }
 
        index_t vertex_first_bundle(index_t v) const {
            return v_first_bndl_[v];
        }
 
        index_t next_around_vertex(index_t bndl) const {
            return bndl_next_around_v_[bndl];
        }
 
        index_t nb_bundles_around_vertex(index_t v) const {
            index_t result = 0;
            for(
                index_t bndl = vertex_first_bundle(v);
                bndl != NO_INDEX;
                bndl = next_around_vertex(bndl)
            ) {
                ++result;
            }
            return result;
        }
 
        index_t opposite(index_t bndl) {
            geo_debug_assert(bndl < nb());
            return (bndl >= nb()/2) ? (bndl-nb()/2) : (bndl+nb()/2);
        }
 
        index_t prev_along_polyline(index_t bndl) {
            index_t v = vertex(bndl,0);
            if(nb_bundles_around_vertex(v) != 2) {
                return NO_INDEX;
            }
            for(
                index_t bndl2 = vertex_first_bundle(v);
                bndl2 != NO_INDEX; bndl2 = next_around_vertex(bndl2)
            ) {
                if(bndl2 != bndl) {
                    return opposite(bndl2);
                }
            }
            geo_assert_not_reached;
        }
 
        index_t next_along_polyline(index_t bndl) {
            index_t v = vertex(bndl,1);
            if(nb_bundles_around_vertex(v) != 2) {
                return NO_INDEX;
            }
            for(
                index_t bndl2 = vertex_first_bundle(v);
                bndl2 != NO_INDEX; bndl2 = next_around_vertex(bndl2)
            ) {
                if(opposite(bndl2) != bndl) {
                    return bndl2;
                }
            }
            geo_assert_not_reached;
        }
 
        bool radial_sort(index_t bndl, RadialSort& RS) {
            geo_debug_assert(bndl < nb());
            if(nb_halfedges(bndl) <= 2) {
                return true;
            }
            auto b = H_.begin() + std::ptrdiff_t(bndl_start_[bndl]);
            auto e = H_.begin() + std::ptrdiff_t(bndl_start_[bndl+1]);
            RS.init(*b);
            std::sort(
                b, e, [&RS](index_t h1, index_t h2) {
                    return RS(h1,h2);
                }
            );
            bool OK = !RS.degenerate();
            bndl_is_sorted_[bndl] = OK;
            return OK;
        }
 
        void set_sorted_halfedges(
            index_t bndl, const vector<index_t>& halfedges
        ) {
            geo_debug_assert(halfedges.size() == nb_halfedges(bndl));
            for(index_t i=0; i<halfedges.size(); ++i) {
                set_halfedge(bndl, i, halfedges[i]);
            }
            bndl_is_sorted_[bndl] = true;
        }
 
        typedef std::pair<index_t, index_t> ChartPos;
 
        void get_sorted_incident_charts(
            index_t bndl, vector<ChartPos>& chart_pos
        );
 
        bool is_sorted(index_t bndl) const {
            geo_debug_assert(bndl < nb());
            return bndl_is_sorted_[bndl];
        }
 
        // private:
        MeshSurfaceIntersection& I_;
        Mesh& mesh_;
        Attribute<index_t> facet_chart_;
        vector<index_t> H_;
        vector<index_t> bndl_start_;
        vector<index_t> v_first_bndl_;
        vector<index_t> bndl_next_around_v_;
        vector<char> bndl_is_sorted_; // not vector<bool>, multithread! (#308)
    } radial_bundles_;
 
    /***************************************************/
 
    class RadialPolylines {
    public:
        RadialPolylines(MeshSurfaceIntersection& I) : I_(I), mesh_(I.mesh_) {
        }
 
        void initialize();
 
        void radial_sort();
 
        index_t nb() const {
            return polyline_start_.size() - 1;
        }
 
        index_as_iterator begin() const {
            return index_as_iterator(0);
        }
 
        index_as_iterator end() const {
            return index_as_iterator(nb());
        }
 
        const_index_ptr_range bundles(index_t polyline) const {
            geo_debug_assert(polyline < nb());
            return const_index_ptr_range(
                B_, polyline_start_[polyline], polyline_start_[polyline+1]
            );
        }
 
        index_t nb_bundles(index_t polyline) const {
            geo_debug_assert(polyline < nb());
            return polyline_start_[polyline+1] - polyline_start_[polyline];
        }
 
        index_t bundle(index_t polyline, index_t li) const {
            geo_debug_assert(polyline < nb());
            geo_debug_assert(li < nb_bundles(polyline));
            return B_[polyline_start_[polyline] + li];
        }
 
        void get_skeleton(Mesh& to, bool trim_fins=false);
 
    private:
        MeshSurfaceIntersection& I_;
        Mesh& mesh_;
        vector<index_t> B_;
        vector<index_t> polyline_start_;
    } radial_polylines_;
    };
 
    /********************************************************************/
 
    enum MeshBooleanOperationFlags {
        MESH_BOOL_OPS_DEFAULT            = 0,
        MESH_BOOL_OPS_VERBOSE            = 1,
        MESH_BOOL_OPS_ATTRIBS            = 2,
        MESH_BOOL_OPS_NO_SIMPLIFY        = 4,
        MESH_BOOL_OPS_NO_CHECK_NEIGHBORS = 8
    };
 
    void GEOGRAM_API mesh_boolean_operation(
        Mesh& result, const Mesh& A, const Mesh& B, const std::string& operation,
        MeshBooleanOperationFlags flags = MESH_BOOL_OPS_DEFAULT
    );
 
    inline void mesh_boolean_operation(
        Mesh& result, const Mesh& A, const Mesh& B, const std::string& operation,
        bool verbose
    ) {
        mesh_boolean_operation(
            result, A, B, operation,
            verbose ? MESH_BOOL_OPS_VERBOSE : MESH_BOOL_OPS_DEFAULT
        );
    }
 
    inline void mesh_union(
        Mesh& result, const Mesh& A, const Mesh& B,
        MeshBooleanOperationFlags flags=MESH_BOOL_OPS_DEFAULT
    ) {
        mesh_boolean_operation(result, A, B, "A+B", flags);
    }
 
    inline void mesh_union(
        Mesh& result, const Mesh& A, const Mesh& B, bool verbose
    ) {
        mesh_boolean_operation(result, A, B, "A+B", verbose);
    }
 
 
    inline void mesh_intersection(
        Mesh& result, const Mesh& A, const Mesh& B,
        MeshBooleanOperationFlags flags=MESH_BOOL_OPS_DEFAULT
    ) {
        mesh_boolean_operation(result, A, B, "A*B", flags);
    }
 
    inline void mesh_intersection(
        Mesh& result, const Mesh& A, const Mesh& B, bool verbose
    ) {
        mesh_boolean_operation(result, A, B, "A*B", verbose);
    }
 
    inline void mesh_difference(
        Mesh& result, const Mesh& A, const Mesh& B,
        MeshBooleanOperationFlags flags=MESH_BOOL_OPS_DEFAULT
    ) {
        mesh_boolean_operation(result, A, B, "A-B", flags);
    }
 
    inline void mesh_difference(
        Mesh& result, const Mesh& A, const Mesh& B, bool verbose
    ) {
        mesh_boolean_operation(result, A, B, "A-B", verbose);
    }
 
    void GEOGRAM_API mesh_remove_intersections(
        Mesh& M, index_t max_iter = 3, bool verbose=false
    );
 
    bool GEOGRAM_API mesh_facets_have_intersection(
        Mesh& M, index_t f1, index_t f2
    );
 
    /**************************************************************************/
}
 
#endif