generic_RVD.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
#define GEOGRAM_VORONOI_GENERIC_RVD
 
#include <geogram/basic/common.h>
#include <geogram/basic/numeric.h>
#include <geogram/voronoi/generic_RVD_utils.h>
#include <geogram/voronoi/RVD_callback.h>
#include <geogram/numerics/predicates.h>
#include <geogram/mesh/index.h>
#include <geogram/basic/geometry_nd.h>
#include <geogram/basic/process.h>
#include <geogram/basic/attributes.h>
#include <geogram/basic/argused.h>
 
#include <deque>
#include <algorithm>
#include <iostream>
 
namespace GEOGen {
 
    template <index_t DIM>
    class RestrictedVoronoiDiagram {
 
        typedef RestrictedVoronoiDiagram<DIM> thisclass;
 
    public:
        static coord_index_t dimension() {
            return DIM;
        }
 
        typedef GEOGen::PointAllocator PointAllocator;
 
        typedef GEOGen::Vertex Vertex;
 
        typedef GEOGen::Polygon Polygon;
 
        typedef GEOGen::ConvexCell Polyhedron;
 
        /********************************************************************/
 
        RestrictedVoronoiDiagram(
            Delaunay* delaunay,
            GEO::Mesh* mesh
        ) :
            mesh_(mesh),
            delaunay_(delaunay),
            intersections_(DIM),
            symbolic_(false),
            check_SR_(true),
            exact_(false)
            {
                delaunay_nn_ = dynamic_cast<GEO::Delaunay_NearestNeighbors*>(
                    delaunay_
                );
                dimension_ = DIM;
                facets_begin_ = UNSPECIFIED_RANGE;
                facets_end_ = UNSPECIFIED_RANGE;
                tets_begin_ = UNSPECIFIED_RANGE;
                tets_end_ = UNSPECIFIED_RANGE;
                connected_components_priority_ = false;
                facet_seed_marking_ = nullptr;
                connected_component_changed_ = false;
                current_connected_component_ = 0;
                cur_stamp_ = NO_INDEX;
                current_facet_ = GEO::max_index_t();
                current_seed_ = GEO::max_index_t();
                current_polygon_ = nullptr;
                current_tet_ = GEO::max_index_t();
                current_polyhedron_ = nullptr;
            }
 
        void set_connected_components_priority(bool x) {
            connected_components_priority_ = x;
        }
 
 
        bool connected_components_priority() const {
            return connected_components_priority_;
        }
 
        const GEO::Mesh* mesh() const {
            return mesh_;
        }
 
        GEO::Mesh* mesh() {
            return mesh_;
        }
 
        const Delaunay* delaunay() const {
            return delaunay_;
        }
 
        Delaunay* delaunay() {
            return delaunay_;
        }
 
        void set_delaunay(Delaunay* delaunay) {
            delaunay_ = delaunay;
            delaunay_nn_ = dynamic_cast<GEO::Delaunay_NearestNeighbors*>(
                delaunay_
            );
        }
 
        void set_mesh(GEO::Mesh* mesh) {
            mesh_ = mesh;
        }
 
        void set_facets_range(index_t facets_begin, index_t facets_end) {
            geo_debug_assert(facets_end >= facets_begin);
            facets_begin_ = facets_begin;
            facets_end_ = facets_end;
        }
 
        void set_tetrahedra_range(index_t tets_begin, index_t tets_end) {
            geo_debug_assert(tets_end >= tets_begin);
            tets_begin_ = tets_begin;
            tets_end_ = tets_end;
        }
 
        index_t nb_facets_in_range() const {
            return facets_end_ - facets_begin_;
        }
 
        index_t nb_tetrahedra_in_range() const {
            return tets_end_ - tets_begin_;
        }
 
        index_t current_facet() const {
            return current_facet_;
        }
 
        index_t current_seed() const {
            return current_seed_;
        }
 
        const Polygon& current_polygon() const {
            return *current_polygon_;
        }
 
        index_t current_tet() const {
            return current_tet_;
        }
 
        const Polyhedron& current_polyhedron() const {
            return *current_polyhedron_;
        }
 
        void set_symbolic(bool x) {
            symbolic_ = x;
            // exact mode requires symbolic mode.
            if(exact_) {
                symbolic_ = true;
            }
        }
 
        bool symbolic() const {
            return symbolic_;
        }
 
        void set_exact_predicates(bool x) {
            exact_ = x;
            // exact mode requires symbolic mode.
            if(exact_) {
                symbolic_ = true;
            }
        }
 
        bool exact_predicates() const {
            return exact_;
        }
 
        void set_check_SR(bool x) {
            check_SR_ = x;
        }
 
        bool check_SR() const {
            return check_SR_;
        }
 
        PointAllocator* point_allocator() {
            return &intersections_;
        }
 
    protected:
        template <class ACTION>
        class PolygonAction {
        public:
            PolygonAction(const ACTION& do_it) :
                do_it_(do_it) {
            }
 
            void operator() (
                index_t v,
                index_t f,
                const Polygon& P
            ) const {
                GEO::geo_argused(f);
                const_cast<ACTION&> ( do_it_)(v, P);
            }
 
        protected:
            const ACTION& do_it_;
        };
 
        template <class ACTION>
        class TriangleAction {
        public:
            TriangleAction(const ACTION& do_it) :
                do_it_(do_it) {
            }
 
            void operator() (
                index_t v,
                index_t f,
                const Polygon& P
            ) const {
                GEO::geo_argused(f);
                for(index_t i = 1; i + 1 < P.nb_vertices(); i++) {
                    const_cast<ACTION&> (do_it_)(
                        v, P.vertex(0), P.vertex(i), P.vertex(i + 1)
                    );
                }
            }
 
        protected:
            const ACTION& do_it_;
        };
 
        template <class ACTION>
        class HalfedgeAction {
        public:
            HalfedgeAction(const ACTION& do_it) :
                do_it_(do_it) {
            }
 
            void operator() (
                index_t v,
                index_t f,
                const Polygon& P
            ) const {
                GEO::geo_argused(v);
                GEO::geo_argused(f);
                for(index_t i = 0; i < P.nb_vertices(); i++) {
                    if(P.vertex(i).check_flag(INTERSECT)) {
                        index_t j = P.next_vertex(i);
                        const_cast<ACTION&> (do_it_)(
                            P.vertex(i), P.vertex(j)
                        );
                    }
                }
            }
 
        protected:
            const ACTION& do_it_;
        };
 
        template <class ACTION>
        class BorderHalfedgeAction {
        public:
            BorderHalfedgeAction(const ACTION& do_it) :
                do_it_(do_it) {
            }
 
            void operator() (
                index_t v,
                index_t f,
                const Polygon& P
            ) const {
                GEO::geo_argused(f);
                for(index_t i = 0; i < P.nb_vertices(); i++) {
                    if(P.vertex(i).check_flag(ORIGINAL)) {
                        if(P.vertex(i).adjacent_facet() == -1) {
                            index_t j = P.next_vertex(i);
                            const_cast<ACTION&> (do_it_)(
                                v, P.vertex(i), P.vertex(j)
                            );
                        }
                    }
                }
            }
 
        private:
            const ACTION& do_it_;
        };
 
        template <class ACTION>
        class PrimalTriangleAction {
        public:
            PrimalTriangleAction(const ACTION& do_it) :
                do_it_(do_it) {
            }
 
            void operator() (
                index_t iv1,
                index_t f,
                const Polygon& P
            ) const {
                GEO::geo_argused(f);
                for(index_t i = 0; i < P.nb_vertices(); i++) {
                    const Vertex& ve = P.vertex(i);
                    // Primal triangles correspond to vertices of
                    // the RVD that are on two bisectors.
                    if(ve.sym().nb_bisectors() == 2) {
                        index_t iv2 = ve.sym().bisector(0);
                        index_t iv3 = ve.sym().bisector(1);
                        // This test generates triangle (iv1,iv2,iv3)
                        // only once (i.e. if iv1 is the vertex with
                        // the smallest index).
                        if(iv1 < iv2 && iv1 < iv3) {
                            const_cast<ACTION&> (do_it_)(iv1, iv2, iv3);
                        }
                    }
                }
            }
 
        protected:
            const ACTION& do_it_;
        };
 
        template <class ACTION>
        class PolyhedronAction {
        public:
            PolyhedronAction(const ACTION& do_it) :
                do_it_(do_it) {
            }
 
            void operator() (
                index_t v,
                index_t t,
                const Polyhedron& C
            ) const {
                const_cast<ACTION&> ( do_it_)(v, t, C);
            }
 
        protected:
            const ACTION& do_it_;
        };
 
        template <class ACTION>
        class VolumetricIntegrationSimplexAction {
        public:
            VolumetricIntegrationSimplexAction(
                const ACTION& do_it,
                bool visit_inner_tets = false,
                bool coherent_triangles = false
            ) :
                do_it_(do_it),
                visit_inner_tets_(visit_inner_tets),
                coherent_triangles_(coherent_triangles)
                {
                }
 
            void operator() (index_t v, index_t t, const Polyhedron& C) const {
                for(index_t cv = 0; cv < C.max_v(); ++cv) {
                    signed_index_t ct = C.vertex_triangle(cv);
                    if(ct == -1) {
                        continue;
                    }
                    geo_debug_assert(C.triangle_is_used(index_t(ct)));
 
                    signed_index_t adjacent = C.vertex_id(cv);
                    signed_index_t v_adj = -1;
                    signed_index_t t_adj = -1;
 
                    if(adjacent < 0) {
                        // Negative adjacent indices correspond to
                        // tet-tet links (ignored when we want to triangulate
                        // the border of the restricted Voronoi cell while
                        // ignoring internal structures).
                        if(!visit_inner_tets_) {
                            continue;
                        }
                        t_adj = -adjacent - 1;
                    } else if(adjacent > 0) {
                        // Positive adjacent indices correspond to
                        // Voronoi seed - Voronoi seed link
                        v_adj = adjacent - 1;
                    }
                    // and adjacent indicex equal to zero corresponds
                    // to tet on border.
 
                    Polyhedron::Corner c1(
                        index_t(ct),
                        index_t(C.find_triangle_vertex(index_t(ct), cv))
                    );
 
                    // If required, ensure that two polygonal facets
                    // seen from two different volumetric cells will
                    // be triangulated coherently.
                    if(coherent_triangles_) {
                        move_to_first_corner_of_facet(C, c1, v);
                    }
 
                    const Vertex& v1 = C.triangle_dual(c1.t);
 
                    Polyhedron::Corner c2 = c1;
                    C.move_to_next_around_vertex(c2);
                    geo_debug_assert(c2 != c1);
 
                    Polyhedron::Corner c3 = c2;
                    C.move_to_next_around_vertex(c3);
                    geo_debug_assert(c3 != c1);
                    do {
                        const Vertex& v2 = C.triangle_dual(c2.t);
                        const Vertex& v3 = C.triangle_dual(c3.t);
                        const_cast<ACTION&> (do_it_)(
                            v, index_t(v_adj), t, index_t(t_adj), v1, v2, v3
                        );
                        c2 = c3;
                        C.move_to_next_around_vertex(c3);
                    } while(c3 != c1);
                }
            }
 
            void move_to_first_corner_of_facet(
                const Polyhedron& C, Polyhedron::Corner& c,
                index_t center_vertex_id
            ) const {
                Polyhedron::Corner first = c;
                Polyhedron::Corner cur = c;
                do {
                    if(symbolic_compare(
                           C.triangle_dual(cur.t),
                           C.triangle_dual(c.t),
                           center_vertex_id
                       )) {
                        c = cur;
                    }
                    C.move_to_next_around_vertex(cur);
                } while(cur != first);
            }
 
            static bool symbolic_compare(
                const Vertex& p1, const Vertex& p2, index_t center_vertex_id
            ) {
                GEO::signed_quadindex K1(
                    signed_index_t(center_vertex_id),
                    p1.sym()[0], p1.sym()[1], p1.sym()[2]
                );
                GEO::signed_quadindex K2(
                    signed_index_t(center_vertex_id),
                    p2.sym()[0], p2.sym()[1], p2.sym()[2]
                );
                return K1 < K2;
            }
 
        protected:
            const ACTION& do_it_;
            bool visit_inner_tets_;
            bool coherent_triangles_;
        };
 
        template <class ACTION>
        class TetrahedronAction {
        public:
            TetrahedronAction(
                const ACTION& do_it
            ) :
                do_it_(do_it)
                {
                }
 
            void operator() (
                index_t v,
                index_t t,
                const Polyhedron& C
            ) const {
 
                // Find a vertex of the current cell,
                // that will be used as the 'origin'
                // vertex
                const Vertex* v0 = nullptr;
                index_t t0;
                for(t0 = 0; t0 < C.max_t(); ++t0) {
                    if(C.triangle_is_used(t0)) {
                        v0 = &C.triangle_dual(t0);
                        break;
                    }
                }
 
                // If current cell is empty, return
                if(v0 == nullptr) {
                    return;
                }
 
                for(index_t cv = 0; cv < C.max_v(); ++cv) {
                    signed_index_t ct = C.vertex_triangle(cv);
                    if(ct == -1) {
                        continue;
                    }
                    geo_debug_assert(C.triangle_is_used(index_t(ct)));
 
                    signed_index_t adjacent = C.vertex_id(cv);
                    signed_index_t v_adj = -1;
                    signed_index_t t_adj = -1;
 
                    if(adjacent < 0) {
                        // Negative adjacent indices correspond to
                        // tet-tet links
                        t_adj = -adjacent - 1;
                    } else if(adjacent > 0) {
                        // Positive adjacent indices correspond to
                        // Voronoi seed - Voroni seed link
                        v_adj = adjacent - 1;
                    }
                    // and adjacent indicex equal to zero corresponds
                    // to tet on border.
 
                    Polyhedron::Corner c1(
                        index_t(ct), C.find_triangle_vertex(index_t(ct), cv)
                    );
 
                    // If the current facet is incident to
                    // the origin vertex, then skip it (else
                    // it would generate flat tetrahedra)
                    if(facet_is_incident_to_vertex(C, c1, t0)) {
                        continue;
                    }
 
                    const Vertex& v1 = C.triangle_dual(c1.t);
 
                    Polyhedron::Corner c2 = c1;
                    C.move_to_next_around_vertex(c2);
                    geo_debug_assert(c2 != c1);
 
                    Polyhedron::Corner c3 = c2;
                    C.move_to_next_around_vertex(c3);
                    geo_debug_assert(c3 != c1);
                    do {
                        const Vertex& v2 = C.triangle_dual(c2.t);
                        const Vertex& v3 = C.triangle_dual(c3.t);
                        const_cast<ACTION&> (do_it_)(
                            v, index_t(v_adj), t, index_t(t_adj),
                            *v0, v1, v2, v3
                        );
                        c2 = c3;
                        C.move_to_next_around_vertex(c3);
                    } while(c3 != c1);
                }
            }
 
        protected:
            bool facet_is_incident_to_vertex(
                const Polyhedron& C, Polyhedron::Corner& c, index_t t
            ) const {
                Polyhedron::Corner first = c;
                Polyhedron::Corner cur = c;
                do {
                    if(cur.t == t) {
                        return true;
                    }
                    C.move_to_next_around_vertex(cur);
                } while(cur != first);
                return false;
            }
 
        protected:
            const ACTION& do_it_;
        };
 
        template <class ACTION>
        class PrimalTetrahedronAction {
        public:
            PrimalTetrahedronAction(const ACTION& do_it) :
                do_it_(do_it) {
            }
 
            void operator() (
                index_t v,
                index_t t,
                const Polyhedron& C
            ) const {
                GEO::geo_argused(t);
                for(index_t it = 0; it < C.max_t(); ++it) {
                    if(C.triangle_is_used(it)) {
                        const SymbolicVertex& sym = C.triangle_dual(it).sym();
                        if(sym.nb_bisectors() == 3) {
                            index_t v1 = sym.bisector(0);
                            index_t v2 = sym.bisector(1);
                            index_t v3 = sym.bisector(2);
                            // This test ensures that the tet (v,v1,v2,v3)
                            // is generated only once.
                            if(v < v1 && v < v2 && v < v3) {
                                const_cast<ACTION&> (do_it_)(v, v1, v2, v3);
                            }
                        }
                    }
                }
            }
 
        protected:
            const ACTION& do_it_;
        };
 
    public:
        template <class ACTION>
        inline void for_each_polygon(const ACTION& action) {
            this->template compute_surfacic<PolygonAction<ACTION> >(
                PolygonAction<ACTION>(action)
            );
        }
 
        template <class TRIACTION>
        inline void for_each_triangle(const TRIACTION& action) {
            this->template compute_surfacic<TriangleAction<TRIACTION> >(
                TriangleAction<TRIACTION>(action)
            );
        }
 
        template <class HEACTION>
        inline void for_each_halfedge(const HEACTION& action) {
            this->template compute_surfacic<HalfedgeAction<HEACTION> >(
                HalfedgeAction<HEACTION>(action)
            );
        }
 
        template <class BOACTION>
        inline void for_each_border_halfedge(const BOACTION& action) {
            this->template compute_surfacic<BorderHalfedgeAction<BOACTION> >(
                BorderHalfedgeAction<BOACTION>(action)
            );
        }
 
        template <class PRIMTRIACTION>
        inline void for_each_primal_triangle(const PRIMTRIACTION& action) {
            bool sym_backup = symbolic();
            set_symbolic(true);
            this->template compute_surfacic<PrimalTriangleAction<PRIMTRIACTION>>(
                PrimalTriangleAction<PRIMTRIACTION>(action)
            );
            set_symbolic(sym_backup);
        }
 
        template <class ACTION>
        inline void for_each_polyhedron(const ACTION& action) {
            this->template compute_volumetric<PolyhedronAction<ACTION> >(
                PolyhedronAction<ACTION>(action)
            );
        }
 
        template <class ACTION>
        inline void for_each_volumetric_integration_simplex(
            const ACTION& action,
            bool visit_inner_tets = false, bool coherent_triangles = false
        ) {
            this->template compute_volumetric<
                VolumetricIntegrationSimplexAction<ACTION>
                >(
                    VolumetricIntegrationSimplexAction<ACTION>(
                        action, visit_inner_tets, coherent_triangles
                    )
                );
        }
 
        template <class ACTION>
        inline void for_each_tetrahedron(
            const ACTION& action
        ) {
            this->template compute_volumetric<TetrahedronAction<ACTION> >(
                TetrahedronAction<ACTION>(
                    action
                )
            );
        }
 
        template <class ACTION>
        inline void for_each_primal_tetrahedron(const ACTION& action) {
            bool sym_backup = symbolic();
            set_symbolic(true);
            this->template compute_volumetric<PrimalTetrahedronAction<ACTION> >(
                PrimalTetrahedronAction<ACTION>(
                    action
                )
            );
            set_symbolic(sym_backup);
        }
 
    protected:
        template <class ACTION>
        inline void compute_surfacic(const ACTION& action) {
            if(connected_components_priority_) {
                this->template compute_surfacic_with_cnx_priority<ACTION>(
                    action
                );
            } else {
                this->template compute_surfacic_with_seeds_priority<ACTION>(
                    action
                );
            }
        }
 
        template <class ACTION>
        inline void compute_surfacic_with_seeds_priority(const ACTION& action) {
            if(
                facets_begin_ == UNSPECIFIED_RANGE &&
                facets_end_ == UNSPECIFIED_RANGE
            ) {
                facets_begin_ = 0;
                facets_end_ = mesh_->facets.nb();
            }
            current_polygon_ = nullptr;
            GEO::vector<index_t> seed_stamp(
                delaunay_->nb_vertices(), index_t(-1)
            );
            GEO::vector<bool> facet_is_marked(facets_end_-facets_begin_, false);
            init_get_neighbors();
 
            FacetSeedStack adjacent_facets;
            SeedStack adjacent_seeds;
            Polygon F;
            GEO::Attribute<double> vertex_weight;
            vertex_weight.bind_if_is_defined(
                mesh_->vertices.attributes(), "weight"
            );
 
            // The algorithm propagates along both the facet-graph of
            // the surface and the 1-skeleton of the Delaunay triangulation,
            // and computes all the relevant intersections between
            // each Voronoi cell and facet.
            for(index_t f = facets_begin_; f < facets_end_; f++) {
                if(!facet_is_marked[f-facets_begin_]) {
                    // Propagate along the facet-graph.
                    facet_is_marked[f-facets_begin_] = true;
                    adjacent_facets.push(
                        FacetSeed(f, find_seed_near_facet(f))
                    );
                    while(!adjacent_facets.empty()) {
                        current_facet_ = adjacent_facets.top().f;
                        current_seed_ = adjacent_facets.top().seed;
                        adjacent_facets.pop();
 
                        // Copy the current facet from the Mesh into
                        // RestrictedVoronoiDiagram's Polygon data structure
                        // (gathers all the necessary information)
                        F.initialize_from_mesh_facet(
                            mesh_, current_facet_, symbolic_, vertex_weight
                        );
 
                        // Propagate along the Delaunay 1-skeleton
                        // This will traverse all the seeds such that their
                        // Voronoi cell has a non-empty intersection with
                        // the current facet.
                        seed_stamp[current_seed_] = current_facet_;
                        adjacent_seeds.push(current_seed_);
 
                        while(!adjacent_seeds.empty()) {
                            current_seed_ = adjacent_seeds.top();
                            adjacent_seeds.pop();
 
                            current_polygon_ = intersect_cell_facet(
                                current_seed_, F
                            );
 
                            action(
                                current_seed_, current_facet_, current_polygon()
                            );
 
                            // Propagate to adjacent facets and adjacent seeds
                            for(index_t v = 0;
                                v < current_polygon().nb_vertices(); v++
                               ) {
                                const Vertex& ve = current_polygon().vertex(v);
                                signed_index_t neigh_f = ve.adjacent_facet();
                                if(
                                    neigh_f >= signed_index_t(facets_begin_) &&
                                    neigh_f < signed_index_t(facets_end_) &&
                                    neigh_f != signed_index_t(current_facet_)
                                ) {
                                    if(!facet_is_marked[
                                           index_t(neigh_f)-facets_begin_
                                       ]) {
                                        facet_is_marked[
                                            index_t(neigh_f)-facets_begin_
                                        ] = true;
                                        adjacent_facets.push(
                                            FacetSeed(
                                                index_t(neigh_f),
                                                current_seed_
                                            )
                                        );
                                    }
                                }
                                signed_index_t neigh_s = ve.adjacent_seed();
                                if(neigh_s != -1) {
                                    if(
                                        seed_stamp[neigh_s] != current_facet_
                                    ) {
                                        seed_stamp[neigh_s] = current_facet_;
                                        adjacent_seeds.push(index_t(neigh_s));
                                    }
                                }
                            }
                        }
                    }
                }
            }
            current_polygon_ = nullptr;
        }
 
        template <class ACTION>
        inline void compute_volumetric(const ACTION& action) {
            if(connected_components_priority_) {
                this->template compute_volumetric_with_cnx_priority<ACTION>(
                    action
                );
            } else {
                this->template compute_volumetric_with_seeds_priority<ACTION>(
                    action
                );
            }
        }
 
        template <class ACTION>
        inline void compute_volumetric_with_seeds_priority(const ACTION& action){
            if(
                tets_begin_ == UNSPECIFIED_RANGE &&
                tets_end_ == UNSPECIFIED_RANGE
            ) {
                tets_begin_ = 0;
                tets_end_ = mesh_->cells.nb();
            }
 
            geo_assert(tets_begin_ != UNSPECIFIED_RANGE);
            geo_assert(tets_end_ != UNSPECIFIED_RANGE);
 
            GEO::vector<index_t> seed_stamp(
                delaunay_->nb_vertices(), index_t(-1)
            );
            GEO::vector<bool> tet_is_marked(tets_end_-tets_begin_, false);
            init_get_neighbors();
 
            TetSeedStack adjacent_tets;
            SeedStack adjacent_seeds;
            Polyhedron C(dimension());
            GEO::Attribute<double> vertex_weight;
            vertex_weight.bind_if_is_defined(
                mesh_->vertices.attributes(), "weight"
            );
 
            current_polyhedron_ = &C;
            // The algorithm propagates along both the facet-graph of
            // the surface and the 1-skeleton of the Delaunay triangulation,
            // and computes all the relevant intersections between
            // each Voronoi cell and facet.
            for(index_t t = tets_begin_; t < tets_end_; ++t) {
                if(!tet_is_marked[t-tets_begin_]) {
                    // Propagate along the tet-graph.
                    tet_is_marked[t-tets_begin_] = true;
                    adjacent_tets.push(
                        TetSeed(t, find_seed_near_tet(t))
                    );
                    while(!adjacent_tets.empty()) {
                        current_tet_ = adjacent_tets.top().f;
                        current_seed_ = adjacent_tets.top().seed;
                        adjacent_tets.pop();
 
                        // Note: current cell could be looked up here,
                        // (from current_tet_) if we chose to keep it
                        // and copy it right before clipping (I am
                        // not sure that it is worth it, lookup time
                        // will be probably fast enough)
 
                        // Propagate along the Delaunay 1-skeleton
                        // This will traverse all the seeds such that their
                        // Voronoi cell has a non-empty intersection with
                        // the current facet.
                        seed_stamp[current_seed_] = current_tet_;
                        adjacent_seeds.push(current_seed_);
 
                        while(!adjacent_seeds.empty()) {
                            current_seed_ = adjacent_seeds.top();
                            adjacent_seeds.pop();
 
                            C.initialize_from_mesh_tetrahedron(
                                mesh_, current_tet_, symbolic_, vertex_weight
                            );
 
                            intersect_cell_cell(
                                current_seed_, C
                            );
 
                            action(
                                current_seed_, current_tet_,
                                current_polyhedron()
                            );
 
                            // Propagate to adjacent tets and adjacent seeds
                            // Iterate on the vertices of the cell (remember:
                            // the cell is represented in dual form)
                            for(index_t v = 0;
                                v < current_polyhedron().max_v(); ++v
                               ) {
 
                                //  Skip clipping planes that are no longer
                                // connected to a cell facet.
                                if(
                                    current_polyhedron().vertex_triangle(v)
                                    == -1
                                ) {
                                    continue;
                                }
 
                                signed_index_t id =
                                    current_polyhedron().vertex_id(v);
                                if(id > 0) {
                                    // Propagate to adjacent seed
                                    index_t neigh_s = index_t(id - 1);
                                    if(seed_stamp[neigh_s] != current_tet_) {
                                        seed_stamp[neigh_s] = current_tet_;
                                        adjacent_seeds.push(neigh_s);
                                    }
                                } else if(id < 0) {
                                    // id==0 corresponds to facet on boundary
                                    //  (skipped)
                                    // id<0 corresponds to adjacent tet index
 
                                    // Propagate to adjacent tet
                                    signed_index_t neigh_t = -id - 1;
                                    if(
                                        neigh_t >=
                                        signed_index_t(tets_begin_) &&
                                        neigh_t <  signed_index_t(tets_end_) &&
                                        neigh_t != signed_index_t(current_tet_)
                                    ) {
                                        if(!tet_is_marked[
                                               index_t(neigh_t)-tets_begin_
                                           ]) {
                                            tet_is_marked[
                                                index_t(neigh_t)-tets_begin_
                                            ] = true;
                                            adjacent_tets.push(
                                                TetSeed(
                                                    index_t(neigh_t),
                                                    current_seed_
                                                )
                                            );
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
            current_polyhedron_ = nullptr;
        }
 
 
        template <class ACTION>
        inline void compute_volumetric_with_cnx_priority(
            const ACTION& action
        ) {
 
            if(
                tets_begin_ == UNSPECIFIED_RANGE &&
                tets_end_ == UNSPECIFIED_RANGE
            ) {
                tets_begin_ = 0;
                tets_end_ = mesh_->cells.nb();
            }
 
            geo_assert(tets_begin_ != UNSPECIFIED_RANGE);
            geo_assert(tets_end_ != UNSPECIFIED_RANGE);
 
            current_polyhedron_ = nullptr;
            init_get_neighbors();
 
            std::deque<TetSeed> adjacent_seeds;
            std::stack<index_t> adjacent_tets;
 
            static constexpr index_t NO_STAMP = index_t(-1);
            GEO::vector<index_t> tet_stamp(
                tets_end_ - tets_begin_, NO_STAMP
            );
 
            TetSeedMarking visited(
                tets_end_ - tets_begin_, delaunay_->nb_vertices()
            );
 
            // Yes, facet_seed_marking_ points to the TetSeedMarking,
            // (TetSeedMarking is typedef-ed as FacetSeedMarking),
            // ugly I know... to be revised.
            facet_seed_marking_ = &visited;
            Polyhedron C(dimension());
            current_polyhedron_ = &C;
 
            current_connected_component_ = 0;
            // index_t C_index = tets_end_ + 1; // Unused (see comment later)
 
            GEO::Attribute<double> vertex_weight;
            vertex_weight.bind_if_is_defined(
                mesh_->vertices.attributes(),"weight"
            );
 
            // The algorithm propagates along both the facet-graph of
            // the surface and the 1-skeleton of the Delaunay triangulation,
            // and computes all the relevant intersections between
            // each Voronoi cell and facet.
            for(index_t t = tets_begin_; t < tets_end_; ++t) {
                if(tet_stamp[t - tets_begin_] == NO_STAMP) {
                    current_tet_ = t;
                    current_seed_ = find_seed_near_tet(t);
 
                    adjacent_seeds.push_back(
                        TetSeed(current_tet_, current_seed_)
                    );
 
                    // Propagate along the Delaunay-graph.
                    while(!adjacent_seeds.empty()) {
                        // Yes, f, because TetSeed is typedef-ed as FacetSeed
                        current_tet_ = adjacent_seeds.front().f;
                        current_seed_ = adjacent_seeds.front().seed;
                        adjacent_seeds.pop_front();
                        if(
                            tet_stamp[current_tet_ - tets_begin_] ==
                            current_seed_
                        ) {
                            continue;
                        }
 
                        if(visited.is_marked(current_tet_, current_seed_)) {
                            continue;
                        }
 
                        connected_component_changed_ = true;
                        adjacent_tets.push(current_tet_);
                        tet_stamp[current_tet_ - tets_begin_] =
                            current_seed_;
                        while(!adjacent_tets.empty()) {
                            current_tet_ = adjacent_tets.top();
                            adjacent_tets.pop();
 
                            // Copy the current tet from the Mesh into
                            // RestrictedVoronoiDiagram's Polyhedron
                            // data structure (gathers all the necessary
                            // information)
 
                            C.initialize_from_mesh_tetrahedron(
                                mesh_, current_tet_, symbolic_, vertex_weight
                            );
 
                            // Note: difference with
                            //    compute_surfacic_with_cnx_priority():
                            // Since intersect_cell_cell() overwrites C, we
                            // need to initialize C from the mesh for each
                            // visited (tet,seed) pair (and the test for
                            // current_tet_ change with C_index is not
                            // used here).
                            // C_index = current_tet_;
 
                            intersect_cell_cell(current_seed_, C);
                            action(
                                current_seed_, current_tet_, current_polyhedron()
                            );
                            connected_component_changed_ = false;
 
                            bool touches_RVC_border = false;
 
                            // Propagate to adjacent tets and adjacent seeds
                            for(index_t v = 0;
                                v < current_polyhedron().max_v(); ++v
                               ) {
 
                                //   Skip clipping planes that are no longer
                                // connected to a cell facet.
                                if(
                                    current_polyhedron().vertex_triangle(v)
                                    == -1
                                ) {
                                    continue;
                                }
 
                                signed_index_t id =
                                    current_polyhedron().vertex_id(v);
                                if(id < 0) {
                                    // id == 0 corresponds to facet on boundary
                                    //   (skipped)
                                    // id < 0 corresponds to adjacent tet index
                                    signed_index_t s_neigh_t = -id-1;
                                    if(
                                        s_neigh_t >= signed_index_t(tets_begin_)
                                        &&
                                        s_neigh_t < signed_index_t(tets_end_)
                                    ) {
                                        geo_debug_assert(
                                            s_neigh_t !=
                                            signed_index_t(current_tet_)
                                        );
                                        index_t neigh_t = index_t(s_neigh_t);
                                        if(
                                            tet_stamp[neigh_t - tets_begin_] !=
                                            current_seed_
                                        ) {
                                            tet_stamp[neigh_t - tets_begin_] =
                                                current_seed_;
                                            adjacent_tets.push(neigh_t);
                                        }
                                    }
                                } else if(id > 0) {
                                    index_t neigh_s = index_t(id-1);
                                    touches_RVC_border = true;
                                    TetSeed ts(current_tet_, neigh_s);
                                    if(!visited.is_marked(ts)) {
                                        adjacent_seeds.push_back(ts);
                                    }
                                }
 
                            }
                            if(touches_RVC_border) {
                                visited.mark(
                                    TetSeed(current_tet_, current_seed_),
                                    current_connected_component_
                                );
                            }
                        }
                        ++current_connected_component_;
                    }
                }
            }
            facet_seed_marking_ = nullptr;
        }
 
 
    public:
        bool facet_seed_is_visited(index_t f, index_t s) const {
            geo_debug_assert(facet_seed_marking_ != nullptr);
            return facet_seed_marking_->is_marked(FacetSeed(f, s));
        }
 
        index_t get_facet_seed_connected_component(index_t f, index_t s) const {
            geo_debug_assert(facet_seed_marking_ != nullptr);
            return facet_seed_marking_->get_connected_component(
                FacetSeed(f, s)
            );
        }
 
        bool connected_component_changed() const {
            return connected_component_changed_;
        }
 
        index_t current_connected_component() const {
            return current_connected_component_;
        }
 
    protected:
        template <class ACTION>
        inline void compute_surfacic_with_cnx_priority(
            const ACTION& action
        ) {
 
            if(
                facets_begin_ == UNSPECIFIED_RANGE &&
                facets_end_ == UNSPECIFIED_RANGE
            ) {
                facets_begin_ = 0;
                facets_end_ = mesh_->facets.nb();
            }
 
            current_polygon_ = nullptr;
            init_get_neighbors();
 
            std::deque<FacetSeed> adjacent_seeds;
            std::stack<index_t> adjacent_facets;
 
            static constexpr index_t NO_STAMP = index_t(-1);
            GEO::vector<index_t> facet_stamp(
                facets_end_ - facets_begin_, NO_STAMP
            );
 
            FacetSeedMarking visited(
                facets_end_ - facets_begin_, delaunay_->nb_vertices()
            );
 
            facet_seed_marking_ = &visited;
            Polygon F;
            current_connected_component_ = 0;
            index_t F_index = facets_end_ + 1;
 
            GEO::Attribute<double> vertex_weight;
            vertex_weight.bind_if_is_defined(
                mesh_->vertices.attributes(),"weight"
            );
 
            // The algorithm propagates along both the facet-graph of
            // the surface and the 1-skeleton of the Delaunay triangulation,
            // and computes all the relevant intersections between
            // each Voronoi cell and facet.
            for(index_t f = facets_begin_; f < facets_end_; ++f) {
                if(facet_stamp[f - facets_begin_] == NO_STAMP) {
                    current_facet_ = f;
                    current_seed_ = find_seed_near_facet(f);
 
                    adjacent_seeds.push_back(
                        FacetSeed(current_facet_, current_seed_)
                    );
 
                    // Propagate along the Delaunay-graph.
                    while(
                        !adjacent_seeds.empty()
                    ) {
                        current_facet_ = adjacent_seeds.front().f;
                        current_seed_ = adjacent_seeds.front().seed;
                        adjacent_seeds.pop_front();
                        if(
                            facet_stamp[current_facet_ - facets_begin_] ==
                            current_seed_
                        ) {
                            continue;
                        }
 
                        if(visited.is_marked(current_facet_, current_seed_)) {
                            continue;
                        }
 
                        connected_component_changed_ = true;
                        adjacent_facets.push(current_facet_);
                        facet_stamp[current_facet_ - facets_begin_] =
                            current_seed_;
                        while(!adjacent_facets.empty()) {
                            current_facet_ = adjacent_facets.top();
                            adjacent_facets.pop();
 
                            // Copy the current facet from the Mesh into
                            // RestrictedVoronoiDiagram's Polygon data structure
                            // (gathers all the necessary information)
                            if(F_index != current_facet_) {
                                F.initialize_from_mesh_facet(
                                    mesh_, current_facet_, symbolic_,
                                    vertex_weight
                                );
                                F_index = current_facet_;
                            }
 
                            current_polygon_ = intersect_cell_facet(
                                current_seed_, F
                            );
                            action(
                                current_seed_, current_facet_, current_polygon()
                            );
                            connected_component_changed_ = false;
 
                            bool touches_RVC_border = false;
 
                            // Propagate to adjacent facets and adjacent seeds
                            for(index_t v = 0;
                                v < current_polygon().nb_vertices(); v++
                               ) {
                                const Vertex& ve = current_polygon().vertex(v);
                                signed_index_t s_neigh_f = ve.adjacent_facet();
                                if(
                                    s_neigh_f >= signed_index_t(facets_begin_)
                                    &&
                                    s_neigh_f < signed_index_t(facets_end_)
                                ) {
                                    geo_debug_assert(
                                        s_neigh_f !=
                                        signed_index_t(current_facet_)
                                    );
                                    index_t neigh_f = index_t(s_neigh_f);
                                    if(
                                        facet_stamp[neigh_f - facets_begin_] !=
                                        current_seed_
                                    ) {
                                        facet_stamp[neigh_f - facets_begin_] =
                                            current_seed_;
                                        adjacent_facets.push(neigh_f);
                                    }
                                }
                                signed_index_t neigh_s = ve.adjacent_seed();
                                if(neigh_s != -1) {
                                    touches_RVC_border = true;
                                    FacetSeed fs(
                                        current_facet_, index_t(neigh_s)
                                    );
                                    if(!visited.is_marked(fs)) {
                                        adjacent_seeds.push_back(fs);
                                    }
                                }
                            }
                            if(touches_RVC_border) {
                                visited.mark(
                                    FacetSeed(current_facet_, current_seed_),
                                    current_connected_component_
                                );
                            }
                        }
                        ++current_connected_component_;
                    }
                }
            }
            facet_seed_marking_ = nullptr;
        }
 
        index_t find_seed_near_facet(index_t f) {
            const double* p = mesh_->vertices.point_ptr(
                mesh_->facets.vertex(f,0)
            );
            return find_seed_near_point(p);
        }
 
        index_t find_seed_near_tet(index_t t) {
            index_t v = mesh_->cells.tet_vertex(t, 0);
            const double* p = mesh_->vertices.point_ptr(v);
            return find_seed_near_point(p);
        }
 
        index_t find_seed_near_point(const double* p) {
            // In order to be compatible with the symbolic
            // perturbation, if the nearest neighbor is
            // non-unique, we need to return the one of
            // lowest index (because in case of several seeds
            // at equal distance, the one of lowest index
            // is guaranteed to have the facet in its Voronoi
            // cell from the point of view of symbolic
            // perturbation).
            if(exact_ && delaunay_nn_ != nullptr) {
                // TODO: may need more than 10
                index_t neighbors[10];
                double neighbors_sq_dist[10];
                index_t nb = 10;
                if(delaunay_nn_->nb_vertices() < nb) {
                    nb = delaunay_nn_->nb_vertices();
                }
                delaunay_nn_->nn_search()->get_nearest_neighbors(
                    nb, p, neighbors, neighbors_sq_dist
                );
                index_t nearest = neighbors[0];
                double min_d = neighbors_sq_dist[0];
                for(index_t i = 1; i < nb; ++i) {
                    if(neighbors_sq_dist[i] != min_d) {
                        break;
                    }
                    if(neighbors[i] < nearest) {
                        nearest = neighbors[i];
                    }
                }
                return nearest;
            }
 
            return delaunay_->nearest_vertex(p);
        }
 
        void swap_polygons(Polygon*& ping, Polygon*& pong) {
            if(ping != &P1 && ping != &P2) {
                // First clipping operation, ping points to F
                // (current facet copied)
                ping = &P2;
                pong = &P1;
            } else {
                std::swap(ping, pong);
            }
        }
 
        Polygon* intersect_cell_facet(index_t seed, Polygon& F) {
            intersections_.clear();
 
            // Initialize ping-pong pointers for Sutherland-Hodgman
            // re-entrant clipping and copy current facet into 'ping' buffer.
            Polygon* ping = &F;
            Polygon* pong = &P2;
 
            // Clip current facet by current Voronoi cell (associated with seed)
            if(delaunay_nn_ != nullptr) {
                clip_by_cell_SR(seed, ping, pong);   // "Security Radius" mode.
            } else {
                clip_by_cell(seed, ping, pong);   // Standard mode.
            }
 
            return ping; // Yes, 'ping', and not 'pong'
                         // see comments in clip_by_cell()
        }
 
        void clip_by_cell_SR(index_t i, Polygon*& ping, Polygon*& pong) {
            // 'Security radius' mode.
            // Note: the vertices of the neighborhood are returned in
            // increasing distance to pi. We stop the clippings as soon as the
            // 'security radius' is reached.
            const double* geo_restrict pi = delaunay_->vertex_ptr(i);
            geo_assume_aligned(pi, geo_dim_alignment(DIM));
 
            index_t jj = 0;
            index_t prev_nb_neighbors = 0;
            neighbors_.resize(0);
            while(neighbors_.size() < delaunay_nn_->nb_vertices() - 1) {
 
                delaunay_nn_->get_neighbors(i, neighbors_);
                if(neighbors_.size() == 0) {
                    return;
                }
                if(prev_nb_neighbors == neighbors_.size()) {
                    return;
                }
 
                for(; jj < neighbors_.size(); jj++) {
                    index_t j = neighbors_[jj];
                    double R2 = 0.0;
                    for(index_t k = 0; k < ping->nb_vertices(); k++) {
                        geo_decl_aligned(double dik);
                        const double* geo_restrict pk = ping->vertex(k).point();
                        geo_assume_aligned(pk, geo_dim_alignment(DIM));
                        dik = GEO::Geom::distance2(pi, pk, dimension());
                        R2 = std::max(R2, dik);
                    }
                    geo_decl_aligned(double dij);
                    const double* geo_restrict pj = delaunay_->vertex_ptr(j);
                    geo_assume_aligned(pj, geo_dim_alignment(DIM));
                    dij = GEO::Geom::distance2(pi, pj, dimension());
                    // A little bit more than 4, because when
                    // exact predicates are used, we need to
                    // include tangent bisectors in the computation.
                    if(dij > 4.1 * R2) {
                        return;
                    }
                    clip_by_plane(*ping, *pong, i, j);
                    swap_polygons(ping, pong);
                }
 
                if(!check_SR_) {
                    return;
                }
 
                index_t nb_neighbors = neighbors_.size();
                prev_nb_neighbors = nb_neighbors;
 
                if(nb_neighbors > 8) {
                    nb_neighbors += nb_neighbors / 8;
                } else {
                    nb_neighbors++;
                }
 
                nb_neighbors = std::min(
                    nb_neighbors,
                    delaunay_nn_->nb_vertices() - 1
                );
 
                delaunay_nn_->enlarge_neighborhood(i, nb_neighbors);
            }
        }
 
        void clip_by_cell(index_t i, Polygon*& ping, Polygon*& pong) {
            get_neighbors(i);
            for(index_t jj = 0; jj < neighbors_.size(); jj++) {
                index_t j = neighbors_[jj];
                clip_by_plane(*ping, *pong, i, j);
                swap_polygons(ping, pong);
            }
        }
 
        void clip_by_plane(
            Polygon& ping, Polygon& pong,
            index_t i, index_t j
        ) {
            ping.clip_by_plane<DIM>(
                pong, intersections_, mesh_, delaunay_, i, j, exact_, symbolic_
            );
        }
 
    public:
        void intersect_cell_cell(index_t seed, Polyhedron& C) {
            // Clip current facet by current Voronoi cell (associated with seed)
            if(delaunay_nn_ != nullptr) {
                clip_by_cell_SR(seed, C);   // "Security Radius" mode.
            } else {
                clip_by_cell(seed, C);   // Standard mode.
            }
        }
 
    protected:
        void clip_by_cell_SR(index_t seed, Polyhedron& C) {
 
            // 'Security radius' mode.
            // Note: the vertices of the neighborhood are returned in
            // increasing distance to pi. We stop the clippings as soon as the
            // 'security radius' is reached.
            const double* geo_restrict pi = delaunay_->vertex_ptr(seed);
            geo_assume_aligned(pi, geo_dim_alignment(DIM));
 
            index_t jj = 0;
            index_t prev_nb_neighbors = 0;
            neighbors_.resize(0);
 
            while(neighbors_.size() < delaunay_nn_->nb_vertices() - 1) {
 
                delaunay_nn_->get_neighbors(seed, neighbors_);
                if(neighbors_.size() == 0) {
                    return;
                }
                if(prev_nb_neighbors == neighbors_.size()) {
                    return;
                }
 
                for(; jj < neighbors_.size(); jj++) {
                    index_t j = neighbors_[jj];
                    double R2 = 0.0;
                    for(index_t k = 0; k < C.max_t(); ++k) {
                        if(!C.triangle_is_used(k)) {
                            continue;
                        }
                        geo_decl_aligned(double dik);
                        const double* geo_restrict pk =
                            C.triangle_dual(k).point();
                        geo_assume_aligned(pk, geo_dim_alignment(DIM));
                        dik = GEO::Geom::distance2(pi, pk, dimension());
                        R2 = std::max(R2, dik);
                    }
                    geo_decl_aligned(double dij);
                    const double* geo_restrict pj = delaunay_->vertex_ptr(j);
                    geo_assume_aligned(pj, geo_dim_alignment(DIM));
                    dij = GEO::Geom::distance2(pi, pj, dimension());
                    // A little bit more than 4, because when
                    // exact predicates are used, we need to
                    // include tangent bisectors in the computation.
                    if(dij > 4.1 * R2) {
                        return;
                    }
                    clip_by_plane(C, seed, j);
                }
 
                if(!check_SR_) {
                    return;
                }
 
                index_t nb_neighbors = neighbors_.size();
                prev_nb_neighbors = nb_neighbors;
 
                if(nb_neighbors > 8) {
                    nb_neighbors += nb_neighbors / 8;
                } else {
                    nb_neighbors++;
                }
 
                nb_neighbors = std::min(
                    nb_neighbors,
                    delaunay_nn_->nb_vertices() - 1
                );
                delaunay_nn_->enlarge_neighborhood(seed, nb_neighbors);
            }
        }
 
        void clip_by_cell(index_t seed, Polyhedron& C) {
            get_neighbors(seed);
            // Check whether cell is empty (may happen with
            // power diagrams)
            if(neighbors_.size() == 0) {
                C.clear();
            }
            for(index_t jj = 0; jj < neighbors_.size(); jj++) {
                index_t j = neighbors_[jj];
                clip_by_plane(C, seed, j);
            }
        }
 
        void clip_by_plane(Polyhedron& C, index_t i, index_t j) {
            C.clip_by_plane<DIM>(
                mesh_, delaunay_, i, j, exact_, symbolic_
            );
        }
 
        void init_get_neighbors() {
            // In dimension 3 (and if we do not used the ANN-based algorithm),
            // we can use the faster 'stamp-based' algorithm for finding the
            // neighbors.
            if(delaunay_->dimension() == 3 && delaunay_->nb_cells() != 0) {
                cur_stamp_ = 0;
                stamp_.assign(delaunay_->nb_vertices(), NO_INDEX);
            }
        }
 
        void get_neighbors(index_t v) {
            if(stamp_.size() == 0) {
                // Used in ANN mode and with higher dimensions.
                delaunay_->get_neighbors(v, neighbors_);
            } else {
                // Used in 3D mode with standard Delaunay.
                // The following loop replaces
                //    delaunay_->get_neighbors(v,neighbors_) ;
                // (and makes the overall algorithm 10 to 30% more efficient)
                neighbors_.resize(0);
                index_t t = index_t(delaunay_->vertex_cell(v));
                do {
                    index_t lv = delaunay_->index(t, v);
                    for(index_t lw = 0; lw < delaunay_->cell_size(); lw++) {
                        if(lw != lv) {
                            index_t w = index_t(delaunay_->cell_vertex(t, lw));
                            if(stamp_[w] != cur_stamp_) {
                                stamp_[w] = cur_stamp_;
                                neighbors_.push_back(w);
                            }
                        }
                    }
                    t = index_t(delaunay_->next_around_vertex(t, lv));
                } while(t != index_t(delaunay_->vertex_cell(v)));
                cur_stamp_++;
            }
        }
 
    protected:
        GEO::Mesh* mesh_;
        Delaunay* delaunay_;
        GEO::Delaunay_NearestNeighbors* delaunay_nn_;
 
        PointAllocator intersections_;
        Polygon* current_polygon_;
        Polygon P1, P2;
        GEO::vector<index_t> neighbors_;
        index_t current_facet_;
        index_t current_seed_;
        Polyhedron* current_polyhedron_;
        index_t current_tet_;
 
        // For optimized get_neighbors().
        index_t cur_stamp_;
        GEO::vector<index_t> stamp_;
 
        bool symbolic_;
        bool check_SR_;
        bool exact_;
 
        coord_index_t dimension_;
 
        static constexpr index_t UNSPECIFIED_RANGE = index_t(-1);
 
        index_t facets_begin_;
        index_t facets_end_;
 
        index_t tets_begin_;
        index_t tets_end_;
 
        bool connected_components_priority_;
        FacetSeedMarking* facet_seed_marking_;
        bool connected_component_changed_;
        index_t current_connected_component_;
 
        RestrictedVoronoiDiagram(const thisclass& rhs) = delete;
 
        thisclass& operator= (const thisclass& rhs) = delete;
    };
}
 
namespace GEO {
 
    typedef GEOGen::SymbolicVertex SymbolicVertex;
}
 
#endif