generic_RVD_polygon.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_VORONOI_GENERIC_RVD_POLYGON
#define GEOGRAM_VORONOI_GENERIC_RVD_POLYGON
 
#include <geogram/basic/common.h>
#include <geogram/voronoi/generic_RVD_vertex.h>
#include <geogram/basic/attributes.h>
 
namespace GEOGen {
 
    class Polygon {
    public:
 
        index_t nb_vertices() const {
            return index_t(vertex_.size());
        }
 
        const Vertex& vertex(index_t i) const {
            geo_debug_assert(i < nb_vertices());
            return vertex_[i];
        }
 
        Vertex& vertex(index_t i) {
            geo_debug_assert(i < nb_vertices());
            return vertex_[i];
        }
 
        index_t next_vertex(index_t i) const {
            geo_debug_assert(i < nb_vertices());
            return
                (i == nb_vertices() - 1) ? 0 : (i + 1)
                ;
        }
 
        index_t prev_vertex(index_t i) const {
            geo_debug_assert(i < nb_vertices());
            return (i == 0) ? (nb_vertices() - 1) : (i - 1);
        }
 
        Vertex* add_vertex(const Vertex& v) {
            vertex_.push_back(v);
            return &*(vertex_.rbegin());
        }
 
        void clear() {
            vertex_.resize(0);
        }
 
        void resize(index_t sz) {
            vertex_.resize(sz);
        }
 
        void initialize_from_mesh_facet(
            const Mesh* mesh, index_t f, bool symbolic,
            const GEO::Attribute<double>& vertex_weight
        );
 
        template <index_t DIM>
        void clip_by_plane(
            Polygon& target, PointAllocator& target_intersections,
            const Mesh* mesh, const Delaunay* delaunay,
            index_t i, index_t j,
            bool exact, bool symbolic
        ) {
            if(exact) {
                clip_by_plane_exact<DIM>(
                    target, target_intersections, mesh, delaunay, i, j
                );
            } else {
                clip_by_plane_fast<DIM>(
                    target, target_intersections, delaunay, i, j, symbolic
                );
            }
        }
 
        void copy(const Polygon& rhs) {
            vertex_ = rhs.vertex_;
        }
 
        void swap(Polygon& rhs) {
            vertex_.swap(rhs.vertex_);
        }
 
    protected:
        template <index_t DIM>
        void clip_by_plane_fast(
            Polygon& target, PointAllocator& target_intersections,
            const Delaunay* delaunay, index_t i, index_t j,
            bool symbolic
        ) const {
            target.clear();
            if(nb_vertices() == 0) {
                return;
            }
 
            const double* geo_restrict pi = delaunay->vertex_ptr(i);
            geo_assume_aligned(pi, geo_dim_alignment(DIM));
            const double* geo_restrict pj = delaunay->vertex_ptr(j);
            geo_assume_aligned(pj, geo_dim_alignment(DIM));
 
            // Compute d = n . m, where n is the
            // normal vector of the bisector [pi,pj]
            // and m the middle point of the bisector.
            geo_decl_aligned(double d);
            d = 0;
            for(coord_index_t c = 0; c < DIM; ++c) {
                d += (pi[c] + pj[c]) * (pi[c] - pj[c]);
            }
 
            // The predecessor of the first vertex is the last vertex
            index_t prev_k = nb_vertices() - 1;
            const Vertex* prev_vk = &(vertex(prev_k));
            const double* geo_restrict prev_pk = prev_vk->point();
            geo_assume_aligned(prev_pk, geo_dim_alignment(DIM));
 
            // We compute:
            //    prev_l = prev_vk . n
            geo_decl_aligned(double prev_l);
            prev_l = 0.0;
            for(coord_index_t c = 0; c < DIM; ++c) {
                prev_l += prev_pk[c] * (pi[c] - pj[c]);
            }
 
            // We compute:
            //    side1(pi,pj,q) = sign(2*q.n - n.m) = sign(2*l - d)
            GEO::Sign prev_status = GEO::geo_sgn(2.0 * prev_l - d);
 
            for(index_t k = 0; k < nb_vertices(); k++) {
                const Vertex* vk = &(vertex(k));
                const double* pk = vk->point();
 
                // We compute: l = vk . n
                geo_decl_aligned(double l);
                l = 0.0;
                for(coord_index_t c = 0; c < DIM; ++c) {
                    l += pk[c] * (pi[c] - pj[c]);
                }
 
                // We compute:
                //   side1(pi,pj,q) = sign(2*q.n - n.m) = sign(2*l - d)
                GEO::Sign status = GEO::geo_sgn(2.0 * l - d);
 
                // If status of edge extremities differ,
                // then there is an intersection.
                if(status != prev_status && (prev_status != 0)) {
                    Vertex I;
                    double* Ipoint = target_intersections.new_item();
                    I.set_point(Ipoint);
                    if(symbolic) {
                        if(
                            !I.sym().intersect_symbolic(
                                prev_vk->sym(), vk->sym(), j
                            )
                        ) {
                            // We encountered a problem. As a workaround,
                            // we copy prev_vk into the result.
                            I = *prev_vk;
                        }
                    }
 
                    // Compute lambda1 and lambda2, the
                    // barycentric coordinates of the intersection I
                    // in the segment [prev_vk vk]
                    // Note that d and l (used for the predicates)
                    // are reused here.
                    double denom = 2.0 * (prev_l - l);
                    double lambda1, lambda2;
 
                    // Shit happens ! [Forrest Gump]
                    if(::fabs(denom) < 1e-20) {
                        lambda1 = 0.5;
                        lambda2 = 0.5;
                    } else {
                        lambda1 = (d - 2.0 * l) / denom;
                        // Note: lambda2 is also given
                        // by (2.0*l2-d)/denom
                        // (but 1.0 - lambda1 is a bit
                        //  faster to compute...)
                        lambda2 = 1.0 - lambda1;
                    }
                    // Compute intersection I by weighting
                    // the edge extremities with the barycentric
                    // coordinates lambda1 and lambda2
                    for(coord_index_t c = 0; c < DIM; ++c) {
                        Ipoint[c] =
                            lambda1 * prev_pk[c] +
                            lambda2 * pk[c];
                    }
                    I.set_weight(
                        lambda1 * prev_vk->weight() + lambda2 * vk->weight()
                    );
                    if(status > 0) {
                        I.copy_edge_from(*prev_vk);
                        I.set_adjacent_seed(signed_index_t(j));
                    } else {
                        I.set_flag(INTERSECT);
                        I.set_adjacent_seed(vk->adjacent_seed());
                    }
                    target.add_vertex(I);
                }
                if(status > 0) {
                    target.add_vertex(*vk);
                }
                prev_vk = vk;
                prev_pk = pk;
                prev_status = status;
                prev_k = k;
                prev_l = l;
            }
        }
 
        template <index_t DIM>
        void clip_by_plane_exact(
            Polygon& target, PointAllocator& target_intersections,
            const Mesh* mesh, const Delaunay* delaunay,
            index_t i, index_t j
        ) {
            target.clear();
            if(nb_vertices() == 0) {
                return;
            }
 
            const double* pi = delaunay->vertex_ptr(i);
            const double* pj = delaunay->vertex_ptr(j);
 
            // The predecessor of the first vertex is the last vertex
            index_t prev_k = nb_vertices() - 1;
            const Vertex* prev_vk = &(vertex(prev_k));
            Sign prev_status = side_exact(
                mesh, delaunay, *prev_vk, pi, pj, DIM
            );
 
            for(index_t k = 0; k < nb_vertices(); ++k) {
                const Vertex* vk = &(vertex(k));
                Sign status = side_exact(mesh, delaunay, *vk, pi, pj, DIM);
 
                // If status of edge extremities differ,
                // there is an intersection.
                if(status != prev_status && (prev_status != 0)) {
 
                    Vertex I;
                    if(
                        !I.sym().intersect_symbolic(
                            prev_vk->sym(), vk->sym(), j
                        )
                    ) {
                        // We encountered a problem. As a workaround,
                        // we copy prev_vk into the result.
                        I = *prev_vk;
                        // geo_assert_not_reached ;
                        // not supposed to happen in exact mode
                    }
                    I.intersect_geom<DIM>(
                        target_intersections, *prev_vk, *vk, pi, pj
                    );
                    if(status > 0) {
                        I.copy_edge_from(*prev_vk);
                        I.set_adjacent_seed(signed_index_t(j));
                    } else {
                        I.set_flag(INTERSECT);
                        I.set_adjacent_seed(vk->adjacent_seed());
                    }
                    target.add_vertex(I);
                }
                if(status > 0) {
                    target.add_vertex(*vk);
                }
                prev_vk = vk;
                prev_status = status;
                prev_k = k;
            }
        }
 
        static Sign side_exact(
            const Mesh* mesh, const Delaunay* delaunay,
            const Vertex& q, const double* pi, const double* pj,
            coord_index_t dim
        );
 
    private:
        GEO::vector<Vertex> vertex_;
    };
}
 
#endif