periodic_delaunay_3d.h

/*
 *  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 PERIODIC_DELAUNAY_TRIANGULATION_3D
#define PERIODIC_DELAUNAY_TRIANGULATION_3D
 
#include <geogram/basic/common.h>
#include <geogram/delaunay/delaunay.h>
#include <geogram/delaunay/periodic.h>
#include <geogram/voronoi/convex_cell.h>
#include <geogram/basic/process.h>
#include <geogram/basic/geometry.h>
#include <stack>
 
#include <geogram/delaunay/delaunay_sync.h>
 
namespace GEO {
 
    class PeriodicDelaunay3dThread;
 
 
    class GEOGRAM_API PeriodicDelaunay3d : public Delaunay, public Periodic {
    public:
 
        struct IncidentTetrahedra {
            std::stack<index_t> S;
            vector<index_t> incident_tets_set;
 
            void clear_incident_tets() {
                incident_tets_set.resize(0);
            }
 
            void add_incident_tet(index_t t) {
                incident_tets_set.push_back(t);
            }
 
            bool has_incident_tet(index_t t) const {
                for(index_t i=0; i<incident_tets_set.size(); ++i) {
                    if(incident_tets_set[i] == t) {
                        return true;
                    }
                }
                return false;
            }
 
            vector<index_t>::const_iterator begin() const {
                return incident_tets_set.begin();
            }
 
            vector<index_t>::const_iterator end() const {
                return incident_tets_set.end();
            }
        };
 
        PeriodicDelaunay3d(bool periodic, double period=1.0);
 
        PeriodicDelaunay3d(const vec3& period);
 
        void set_vertices(
            index_t nb_vertices, const double* vertices
        ) override;
 
        void set_weights(const double* weights);
 
        void compute();
 
        void use_exact_predicates_for_convex_cell(bool x) {
            convex_cell_exact_predicates_ = x;
        }
 
        vec3 vertex(index_t v) const {
            if(!periodic_) {
                geo_debug_assert(v < nb_vertices());
                return vec3(vertices_ + 3*v);
            }
            index_t instance = v/nb_vertices_non_periodic_;
            v = v%nb_vertices_non_periodic_;
            vec3 result(vertices_ + 3*v);
            result.x += double(translation[instance][0]) * period_.x;
            result.y += double(translation[instance][1]) * period_.y;
            result.z += double(translation[instance][2]) * period_.z;
            return result;
        }
 
        double weight(index_t v) const {
            if(weights_ == nullptr) {
                return 0.0;
            }
            return periodic_ ? weights_[periodic_vertex_real(v)] : weights_[v] ;
        }
 
        index_t nearest_vertex(const double* p) const override;
 
        void set_BRIO_levels(const vector<index_t>& levels) override;
 
        void get_incident_tets(index_t v, IncidentTetrahedra& W) const;
 
        void copy_Laguerre_cell_from_Delaunay(
            GEO::index_t i,
            ConvexCell& C,
            IncidentTetrahedra& W
        ) const;
 
        void copy_Laguerre_cell_from_Delaunay(
            GEO::index_t i,
            ConvexCell& C
        ) const {
            IncidentTetrahedra W;
            copy_Laguerre_cell_from_Delaunay(i,C,W);
        }
 
        bool has_empty_cells() const {
            return has_empty_cells_;
        }
 
        void save_cells(const std::string& basename, bool clipped);
 
    protected:
 
        GEO::index_t copy_Laguerre_cell_facet_from_Delaunay(
            GEO::index_t i,
            const GEO::vec3& Pi,
            double wi,
            double Pi_len2,
            GEO::index_t t,
            ConvexCell& C,
            IncidentTetrahedra& W
        ) const;
 
 
        index_t compress(bool shrink=true);
 
        void update_v_to_cell() override;
 
        void update_cicl() override;
 
        void handle_periodic_boundaries();
 
 
        void handle_periodic_boundaries_phase_I();
 
        bool Laguerre_vertex_is_in_conflict_with_plane(index_t t, vec4 P) const;
 
        void handle_periodic_boundaries_phase_II();
 
        void insert_vertices(const char* phase, index_t b, index_t e);
 
        void insert_vertices_with_BRIO(
            const char* phase, const vector<index_t>& levels
        );
 
        void check_volume();
 
        PeriodicDelaunay3dThread* thread(index_t t) {
            geo_debug_assert(t < threads_.size());
            return reinterpret_cast<PeriodicDelaunay3dThread*>(
                threads_[t].get()
            );
        }
 
        index_t nb_threads() const {
            return index_t(threads_.size());
        }
 
        void check_max_t();
 
    private:
        friend class PeriodicDelaunay3dThread;
 
        bool periodic_;
        vec3 period_;
 
        const double* weights_;
        vector<signed_index_t> cell_to_v_store_;
        vector<signed_index_t> cell_to_cell_store_;
        vector<index_t> cell_next_;
 
        CellStatusArray cell_status_;
 
        ThreadGroup threads_;
        vector<index_t> reorder_;
        vector<index_t> levels_;
 
        bool debug_mode_;
 
        bool verbose_debug_mode_;
 
        bool benchmark_mode_;
 
 
        vector<Numeric::uint32> vertex_instances_;
 
        bool update_periodic_v_to_cell_;
        vector<index_t> periodic_v_to_cell_rowptr_;
        vector<index_t> periodic_v_to_cell_data_;
 
        bool has_empty_cells_;
 
        index_t nb_reallocations_;
 
        bool convex_cell_exact_predicates_;
    };
 
 
}
 
#endif