mesh_surface_intersection_internal.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_INTERNAL
#define GEOGRAM_MESH_MESH_SURFACE_INTERSECTION_INTERNAL
 
#include <geogram/basic/common.h>
#include <geogram/mesh/mesh_surface_intersection.h>
#include <geogram/mesh/mesh.h>
#include <geogram/mesh/mesh_io.h>
#include <geogram/mesh/index.h>
#include <geogram/mesh/triangle_intersection.h>
#include <geogram/delaunay/CDT_2d.h>
#include <geogram/numerics/exact_geometry.h>
 
#include <map>
 
namespace GEO {
 
    class MeshInTriangle : public CDTBase2d {
    public:
 
        typedef exact::vec3h ExactPoint;
 
        /***************************************************************/
 
        class Edge {
        public:
            Edge(
                index_t v1_in = NO_INDEX,
                index_t v2_in = NO_INDEX,
                index_t f2    = NO_INDEX,
                TriangleRegion R2 = T2_RGN_T
            ) : v1(v1_in),
                v2(v2_in) {
                sym.f2 = f2;
                sym.R2 = R2;
            }
            index_t v1,v2; // The two extremities of the edge
            struct {       // Symbolic information: this edge = f1 /\ f2.R2
                index_t        f2;
                TriangleRegion R2;
            } sym;
        };
 
        /***************************************************************/
 
        class Vertex {
        public:
 
            enum Type {
                UNINITIALIZED, MESH_VERTEX, PRIMARY_ISECT, SECONDARY_ISECT
            };
 
            Vertex(MeshInTriangle* M, index_t f, index_t lv) {
                geo_debug_assert(f == M->f1_);
                type = MESH_VERTEX;
                mit = M;
                init_sym(f, NO_INDEX, TriangleRegion(lv), T2_RGN_T);
                init_geometry(compute_geometry());
            }
 
            Vertex(
                MeshInTriangle* M,
                index_t f1, index_t f2,
                TriangleRegion R1, TriangleRegion R2
            ) {
                geo_debug_assert(f1 == M->f1_);
                type = PRIMARY_ISECT;
                mit = M;
                init_sym(f1,f2,R1,R2);
                init_geometry(compute_geometry());
            }
 
            Vertex(
                MeshInTriangle* M, const ExactPoint& point_exact_in
            ) {
                type = SECONDARY_ISECT;
                mit = M;
                init_sym(NO_INDEX, NO_INDEX, T1_RGN_T, T2_RGN_T);
                init_geometry(point_exact_in);
            }
 
            Vertex() {
                type = UNINITIALIZED;
                mit = nullptr;
                init_sym(NO_INDEX, NO_INDEX, T1_RGN_T, T2_RGN_T);
                mesh_vertex_index = NO_INDEX;
            }
 
            const Mesh& mesh() const {
                return mit->mesh();
            }
 
            void print(std::ostream& out=std::cerr) const;
 
            std::string to_string() const {
                std::ostringstream out;
                print(out);
                return out.str();
            }
 
            vec2 get_UV_approx() const {
                double u = point_exact[mit->u_].estimate();
                double v = point_exact[mit->v_].estimate();
                double w = point_exact.w.estimate();
                return vec2(u/w,v/w);
            }
 
        protected:
 
            void init_sym(
                index_t f1, index_t f2, TriangleRegion R1, TriangleRegion R2
            ) {
                sym.f1 = f1;
                sym.f2 = f2;
                sym.R1 = R1;
                sym.R2 = R2;
                mesh_vertex_index = NO_INDEX;
            }
 
            ExactPoint compute_geometry();
 
            void init_geometry(const ExactPoint& P);
 
        public:
            MeshInTriangle* mit;
            ExactPoint point_exact; // Exact homogeneous coords using expansions
            Type type;          // MESH_VERTEX, PRIMARY_ISECT or SECONDARY_ISECT
            index_t mesh_vertex_index; // Global mesh vertex index once created
            struct {                   // Symbolic information - tri-tri isect
                index_t f1,f2;         //   global facet indices in mesh
                TriangleRegion R1,R2;  //   triangle regions
            } sym;
#ifndef GEOGRAM_USE_EXACT_NT
            double l; // precomputed approximated (p[u]^2 + p[v]^2) / p.w^2
#endif
        };
 
        /***************************************************************/
 
        MeshInTriangle(MeshSurfaceIntersection& EM);
 
        const Mesh& mesh() const {
            return mesh_;
        }
 
        Mesh& target_mesh() {
            return exact_mesh_.target_mesh();
        }
 
        void set_dry_run(bool x) {
            dry_run_ = x;
        }
 
        void save_constraints(const std::string& filename) {
            Mesh M;
            get_constraints(M);
            mesh_save(M,filename);
        }
 
        void begin_facet(index_t f);
 
        index_t add_vertex(index_t f2, TriangleRegion R1, TriangleRegion R2);
 
        void add_edge(
            index_t f2,
            TriangleRegion AR1, TriangleRegion AR2,
            TriangleRegion BR1, TriangleRegion BR2
        );
 
        void commit();
 
 
        void clear() override;
 
    protected:
        void get_constraints(Mesh& M, bool with_edges=true) const;
 
        vec3 mesh_vertex(index_t v) const {
            return vec3(mesh().vertices.point_ptr(v));
        }
 
        vec3 mesh_facet_vertex(index_t f, index_t lv) const {
            index_t v = mesh().facets.vertex(f,lv);
            return mesh_vertex(v);
        }
 
        vec2 mesh_vertex_UV(index_t v) const {
            const double* p = mesh().vertices.point_ptr(v);
            return vec2(p[u_], p[v_]);
        }
 
        vec2 mesh_facet_vertex_UV(index_t f, index_t lv) const {
            index_t v = mesh().facets.vertex(f,lv);
            return mesh_vertex_UV(v);
        }
 
 
        void log_err() const {
            std::cerr << "Houston, we got a problem (while remeshing facet "
                      << f1_ << "):" << std::endl;
        }
 
    protected:
 
        /********************** CDTBase2d overrides ***********************/
 
        Sign orient2d(index_t v1,index_t v2,index_t v3) const override;
 
        Sign incircle(
            index_t v1,index_t v2,index_t v3,index_t v4
        ) const override;
 
        index_t create_intersection(
            index_t e1, index_t i, index_t j,
            index_t e2, index_t k, index_t l
        ) override;
 
        void get_edge_edge_intersection(
            index_t e1, index_t e2, ExactPoint& I
        ) const;
 
        void get_edge_edge_intersection_2D(
            index_t e1, index_t e2, ExactPoint& I
        ) const;
 
    public:
        void save(const std::string& filename) const override;
 
    protected:
        void begin_insert_transaction() override;
        void commit_insert_transaction() override;
        void rollback_insert_transaction() override;
 
    private:
        MeshSurfaceIntersection& exact_mesh_;
        const Mesh& mesh_;
        index_t f1_;
        index_t latest_f2_;
        index_t latest_f2_count_;
        coord_index_t f1_normal_axis_;
        coord_index_t u_; // = (f1_normal_axis_ + 1)%3
        coord_index_t v_; // = (f1_normal_axis_ + 2)%3
        vector<Vertex> vertex_;
        vector<Edge> edges_;
        bool has_planar_isect_;
        bool dry_run_;
        mutable std::map<trindex, Sign> pred_cache_;
        bool use_pred_cache_insert_buffer_;
        mutable std::vector< std::pair<trindex, Sign> >
        pred_cache_insert_buffer_;
    };
 
    /*************************************************************************/
 
    struct IsectInfo {
    public:
 
        void flip() {
            std::swap(f1,f2);
            A_rgn_f1 = swap_T1_T2(A_rgn_f1);
            A_rgn_f2 = swap_T1_T2(A_rgn_f2);
            std::swap(A_rgn_f1, A_rgn_f2);
            B_rgn_f1 = swap_T1_T2(B_rgn_f1);
            B_rgn_f2 = swap_T1_T2(B_rgn_f2);
            std::swap(B_rgn_f1, B_rgn_f2);
        }
 
        bool is_point() const {
            return
                A_rgn_f1 == B_rgn_f1 &&
                A_rgn_f2 == B_rgn_f2 ;
        }
 
        index_t f1;
        index_t f2;
        TriangleRegion A_rgn_f1;
        TriangleRegion A_rgn_f2;
        TriangleRegion B_rgn_f1;
        TriangleRegion B_rgn_f2;
    };
 
    /**********************************************************************/
 
    class CoplanarFacets {
    public:
        static constexpr index_t NON_MANIFOLD = index_t(-2);
        typedef MeshSurfaceIntersection::ExactPoint ExactPoint;
 
        CoplanarFacets(
            MeshSurfaceIntersection& I, bool clear_attributes,
            double angle_tolerance = 0.0
        );
 
        void get(index_t f, index_t group_id);
 
        void mark_vertices_to_keep();
 
        void save_borders(const std::string& filename);
 
        void save_facet_group(const std::string& filename);
 
        void triangulate();
 
    protected:
 
        void find_coplanar_facets();
 
        bool triangles_are_coplanar(
            const vec3& p1, const vec3& p2, const vec3& p3,
            const vec3& q1, const vec3& q2, const vec3& q3
        ) const;
 
 
        bool edges_are_colinear(
            const ExactPoint& P1, const ExactPoint& P2, const ExactPoint& P3
        ) const;
 
 
 
    public:
        ExactCDT2d CDT;
 
        index_t nb_facets() {
            return facets_.size();
        }
 
        void mark_facets(vector<index_t>& facet_is_marked) {
            for(index_t f: facets_) {
                facet_is_marked[f] = 1;
            }
        }
 
    private:
        MeshSurfaceIntersection& I_;
        Mesh& mesh_;
        const Mesh& mesh_copy_;
        double angle_tolerance_;
        index_t group_id_;
        Attribute<index_t> facet_group_;
        Attribute<bool> keep_vertex_;
        Attribute<bool> c_is_coplanar_;
        Attribute<bool> f_is_flipped_;
        vector<bool>    f_visited_;
        vector<bool>    h_visited_;
        vector<bool>    v_visited_;
        vector<index_t> v_idx_;
        coord_index_t   u_;
        coord_index_t   v_;
 
        /***********************************************************/
 
        vector<index_t> vertices_;
        vector<index_t> facets_;
 
        class Halfedges {
        public:
 
            Halfedges(
                CoplanarFacets& coplanar_facets
            ) : mesh_(coplanar_facets.mesh_) {
            }
 
            void initialize() {
                // Use resize rather than assign so that we do not traverse
                // all the halfedges of the mesh_
                v_first_halfedge_.resize(mesh_.vertices.nb(), NO_INDEX);
                h_next_around_v_.resize(mesh_.facet_corners.nb(), NO_INDEX);
                // We only need to reset the halfedges of this set of coplanar
                // facets.
                for(index_t h: halfedges_) {
                    v_first_halfedge_[vertex(h,0)] = NO_INDEX;
                    h_next_around_v_[h] = NO_INDEX;
                }
                halfedges_.resize(0);
            }
 
            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);
            }
 
            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;
            }
 
            void add(index_t h) {
                halfedges_.push_back(h);
                index_t v1 = vertex(h,0);
                h_next_around_v_[h] = v_first_halfedge_[v1];
                v_first_halfedge_[v1] = h;
            }
 
            vector<index_t>::const_iterator begin() const {
                return halfedges_.begin();
            }
 
            vector<index_t>::const_iterator end() const {
                return halfedges_.end();
            }
 
            index_t vertex_first_halfedge(index_t v) const {
                return v_first_halfedge_[v];
            }
 
            index_t next_around_vertex(index_t h) const {
                return h_next_around_v_[h];
            }
 
            index_t nb_halfedges_around_vertex(index_t v) const {
                index_t result = 0;
                for(
                    index_t h = vertex_first_halfedge(v);
                    h != NO_INDEX;
                    h = next_around_vertex(h)
                ) {
                    ++result;
                }
                return result;
            }
 
            index_t next_along_polyline(index_t h) const {
                index_t v2 = vertex(h,1);
                if(nb_halfedges_around_vertex(v2) != 1) {
                    return NO_INDEX;
                }
                return vertex_first_halfedge(v2);
            }
 
        private:
            Mesh& mesh_;
            vector<index_t> halfedges_;
            vector<index_t> v_first_halfedge_;
            vector<index_t> h_next_around_v_;
        } halfedges_;
 
 
        class Polylines {
        public:
 
            Polylines(CoplanarFacets& CF) : CF_(CF) {
            }
 
            void initialize() {
                H_.resize(0);
                polyline_start_.resize(0);
                polyline_start_.push_back(0);
            }
 
            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 halfedges(index_t polyline) const {
                geo_debug_assert(polyline < nb());
                return const_index_ptr_range(
                    H_, polyline_start_[polyline], polyline_start_[polyline+1]
                );
            }
 
            void begin_polyline() {
            }
 
            void end_polyline() {
                polyline_start_.push_back(H_.size());
            }
 
            void add_halfedge(index_t h) {
                H_.push_back(h);
            }
 
            index_t first_vertex(index_t polyline) const {
                index_t h = H_[polyline_start_[polyline]];
                return CF_.halfedges_.vertex(h,0);
            }
 
            index_t last_vertex(index_t polyline) const {
                index_t h = H_[polyline_start_[polyline+1]-1];
                return CF_.halfedges_.vertex(h,1);
            }
 
            index_t prev_first_vertex(index_t polyline) const {
                if(first_vertex(polyline) != last_vertex(polyline)) {
                    return NO_INDEX;
                }
                index_t h = H_[polyline_start_[polyline+1]-1];
                return CF_.halfedges_.vertex(h,0);
            }
 
        private:
            CoplanarFacets& CF_;
            vector<index_t> H_;
            vector<index_t> polyline_start_;
        } polylines_;
 
    };
 
    /**********************************************************************/
 
}
 
#endif