convex_cell.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_CONVEX_CELL
#define GEOGRAM_VORONOI_CONVEX_CELL
 
#ifndef STANDALONE_CONVEX_CELL
#include <geogram/basic/common.h>
#include <geogram/basic/memory.h>
#include <geogram/basic/numeric.h>
#include <geogram/basic/geometry.h>
#  ifndef GEOGRAM_PSM
#  include <geogram/basic/attributes.h>
#  endif
#endif
 
#include <string>
#include <vector>
#include <iostream>
#include <cmath>
#include <cassert>
 
 
 
#ifndef STANDALONE_CONVEX_CELL
namespace GEO {
    class Mesh;
    class PeriodicDelaunay3d;
}
#endif
 
 
namespace VBW {
 
#ifdef STANDALONE_CONVEX_CELL
    using std::vector;
    typedef unsigned int index_t;
    typedef unsigned int global_index_t;
#   define vbw_assert(x) assert(x)
    struct vec2 {
        double x;
        double y;
    };
    struct vec3 {
        double x;
        double y;
        double z;
    };
    struct vec4 {
        double x;
        double y;
        double z;
        double w;
    };
#else
    using GEO::vector;
    typedef unsigned int index_t;        // Always 32 bits
    typedef GEO::index_t global_index_t; // Possibly 64 bits in GARGANTUA mode
#   define vbw_assert(x) geo_debug_assert(x)
    using GEO::vec2;
    using GEO::vec3;
    using GEO::vec4;
#endif
 
/******************************************************************************/
 
 
    inline vec2 make_vec2(
        double x, double y
    ) {
        vec2 result;
        result.x = x;
        result.y = y;
        return result;
    }
 
 
    inline vec3 make_vec3(
        double x, double y, double z
    ) {
        vec3 result;
        result.x = x;
        result.y = y;
        result.z = z;
        return result;
    }
 
 
    inline vec3 cross(vec3 v1, vec3 v2) {
        return make_vec3(
            v1.y*v2.z - v1.z*v2.y,
            v1.z*v2.x - v1.x*v2.z,
            v1.x*v2.y - v1.y*v2.x
        );
    }
 
    inline double dot(vec3 v1, vec3 v2) {
        return (
            v1.x*v2.x + v1.y*v2.y + v1.z*v2.z
        );
    }
 
    inline double squared_length(vec3 v) {
        return (v.x*v.x + v.y*v.y + v.z*v.z);
    }
 
    inline double squared_distance(vec3 v, vec3 w) {
        double dx = w.x-v.x;
        double dy = w.y-v.y;
        double dz = w.z-v.z;
        return (dx*dx+dy*dy+dz*dz);
    }
 
    inline double length(vec3 v) {
        return ::sqrt(squared_length(v));
    }
 
    inline vec3 normalize(vec3 v) {
        double s = 1.0/length(v);
        return make_vec3(
            s*v.x, s*v.y, s*v.z
        );
    }
 
    inline vec4 make_vec4(
        double x, double y, double z, double w
    ) {
        vec4 result;
        result.x = x;
        result.y = y;
        result.z = z;
        result.w = w;
        return result;
    }
 
    inline double dot(vec4 v1, vec4 v2) {
        return (
            v1.x*v2.x + v1.y*v2.y +
            v1.z*v2.z + v1.w*v2.w
        );
    }
 
    inline double squared_length(vec4 v) {
        return (
            v.x*v.x + v.y*v.y +
            v.z*v.z + v.w*v.w
        );
    }
 
    inline double length(vec4 v) {
        return ::sqrt(squared_length(v));
    }
 
    inline double squared_point_plane_distance(VBW::vec3 p, VBW::vec4 P) {
        double result = P.x*p.x + P.y*p.y + P.z*p.z + P.w;
        result = (result*result) / (P.x*P.x + P.y*P.y + P.z*P.z);
        return result;
    }
 
    enum {
        CONFLICT_MASK  = 32768, 
        MARKED_MASK    = 16384, 
        END_OF_LIST    = 16383, 
        VERTEX_AT_INFINITY = 0  
    };
 
 
    typedef unsigned char uchar;
 
    typedef unsigned short ushort;
 
    struct Triangle {
        ushort i;
        ushort j;
        ushort k;
        ushort operator[](unsigned int index) const {
            vbw_assert(index < 3);
            return (&i)[index];
        }
        ushort& operator[](unsigned int index) {
            vbw_assert(index < 3);
            return (&i)[index];
        }
    };
 
    inline Triangle make_triangle(
        ushort i, ushort j, ushort k
    ) {
        Triangle result;
        result.i = i;
        result.j = j;
        result.k = k;
        return result;
    }
 
    struct TriangleWithFlags : public Triangle {
        ushort flags;
    };
 
    inline TriangleWithFlags make_triangle_with_flags(
        ushort i, ushort j, ushort k, ushort f
    ) {
        TriangleWithFlags result;
        result.i = i;
        result.j = j;
        result.k = k;
        result.flags = f;
        return result;
    }
 
 
/******************************************************************************/
 
    inline double det2x2(
        double a11, double a12,
        double a21, double a22
    ) {
        return a11*a22 - a12*a21;
    }
 
    inline double det3x3(
        double a11, double a12, double a13,
        double a21, double a22, double a23,
        double a31, double a32, double a33
    ) {
        return
            a11*det2x2(a22,a23,a32,a33)
            -a21*det2x2(a12,a13,a32,a33)
            +a31*det2x2(a12,a13,a22,a23);
    }
 
    inline double det4x4(
        double a11, double a12, double a13, double a14,
        double a21, double a22, double a23, double a24,
        double a31, double a32, double a33, double a34,
        double a41, double a42, double a43, double a44
    ) {
        double m12 = a21*a12 - a11*a22;
        double m13 = a31*a12 - a11*a32;
        double m14 = a41*a12 - a11*a42;
        double m23 = a31*a22 - a21*a32;
        double m24 = a41*a22 - a21*a42;
        double m34 = a41*a32 - a31*a42;
 
        double m123 = m23*a13 - m13*a23 + m12*a33;
        double m124 = m24*a13 - m14*a23 + m12*a43;
        double m134 = m34*a13 - m14*a33 + m13*a43;
        double m234 = m34*a23 - m24*a33 + m23*a43;
 
        return (m234*a14 - m134*a24 + m124*a34 - m123*a44);
    }
 
/******************************************************************************/
 
    enum ConvexCellFlag {
        None        = 0, 
        WithVGlobal = 1, 
        WithTFlags  = 2  
    };
 
    typedef index_t ConvexCellFlags;
 
    class GEOGRAM_API ConvexCell {
    public:
 
    ConvexCell(ConvexCellFlags flags = None);
 
#ifndef STANDALONE_CONVEX_CELL
    void use_exact_predicates(bool x) {
        use_exact_predicates_ = x;
    }
#endif
 
    bool has_vglobal() const {
        return has_vglobal_;
    }
 
    bool has_tflags() const {
        return has_tflags_;
    }
 
    void create_vglobal() {
        if(!has_vglobal()) {
            has_vglobal_ = true;
            vglobal_.assign(max_v(), global_index_t(-1));
        }
    }
 
    void clear();
 
    void init_with_box(
        double xmin, double ymin, double zmin,
        double xmax, double ymax, double zmax
    );
 
    void init_with_tet(
        vec4 P0, vec4 P1, vec4 P2, vec4 P3
    );
 
    void init_with_tet(
        vec4 P0, vec4 P1, vec4 P2, vec4 P3,
        global_index_t P0_global_index,
        global_index_t P1_global_index,
        global_index_t P2_global_index,
        global_index_t P3_global_index
    );
 
    void save(const std::string& filename, double shrink=0.0) const;
 
 
    index_t save(
        std::ostream& out, global_index_t v_offset=1, double shrink=0.0,
        bool borders_only=false
    ) const;
 
#if !defined(STANDALONE_CONVEX_CELL) && !defined(GEOGRAM_PSM)
    void append_to_mesh(
        GEO::Mesh* mesh,
        double shrink=0.0, bool borders_only=false,
        GEO::Attribute<GEO::index_t>* facet_attr=nullptr
    ) const;
 
#endif
 
    void for_each_Voronoi_vertex(
        index_t v,
        std::function<void(index_t)> vertex
    );
 
    void clip_by_plane(vec4 P);
 
    void clip_by_plane(vec4 P, global_index_t j);
 
 
    void clip_by_plane(
        vec4 P, global_index_t P_global_index,
        std::function<bool(ushort,ushort)> triangle_conflict_predicate
    );
 
    void clip_by_plane_fast(vec4 P);
 
    void clip_by_plane_fast(vec4 P, global_index_t j);
 
    index_t nb_t() const {
        return nb_t_;
    }
 
    index_t nb_v() const {
        return nb_v_;
    }
 
    index_t create_vertex(vec4 P) {
        if(nb_v_ == max_v_) {
            grow_v();
        }
        plane_eqn_[nb_v_] = P;
        index_t result = nb_v_;
        ++nb_v_;
        return result;
    }
 
    index_t create_vertex(vec4 P, global_index_t v) {
        index_t result = create_vertex(P);
        vglobal_[nb_v()-1] = v;
        return result;
    }
 
    index_t create_triangle(index_t i, index_t j, index_t k) {
        vbw_assert(i < nb_v());
        vbw_assert(j < nb_v());
        vbw_assert(k < nb_v());
        return new_triangle(i,j,k);
    }
 
    void kill_vertex(index_t v);
 
    bool vertex_is_contributing(index_t v) const {
        if(!geometry_dirty_) {
            return v2t_[v] != END_OF_LIST;
        }
        index_t t = first_valid_;
        while(t != END_OF_LIST) {
            TriangleWithFlags T = get_triangle_and_flags(t);
            if(T.i == v || T.j == v || T.k == v) {
                return true;
            }
            t = index_t(T.flags);
        }
        return false;
    }
 
    index_t vertex_triangle(index_t v) const {
        geo_assert(!geometry_dirty_);
        return v2t_[v];
    }
 
    void compute_geometry();
 
    double facet_area(index_t v) const;
 
    double volume() const;
 
    vec3 barycenter() const;
 
    void compute_mg(double& m, vec3& mg) const ;
 
 
    double squared_radius(vec3 center) const;
 
    double squared_inner_radius(vec3 center) const;
 
 
    bool empty() const {
        return first_valid_ == END_OF_LIST;
    }
 
    global_index_t v_global_index(index_t lv) const {
        vbw_assert(has_vglobal_);
        vbw_assert(lv < nb_v());
        return vglobal_[lv];
    }
 
    void set_v_global_index(index_t lv, global_index_t v) {
        vbw_assert(has_vglobal_);
        vbw_assert(lv < nb_v());
        vglobal_[lv] = v;
    }
 
    bool has_v_global_index(global_index_t v) const;
 
    ushort first_triangle() const {
        return ushort(first_valid_);
    }
 
    ushort next_triangle(ushort t) const {
        return get_triangle_flags(t);
    }
 
    vec3 triangle_point(ushort t) const {
        if(geometry_dirty_) {
            vec4 result = compute_triangle_point(t);
            vbw_assert(result.w != 0.0);
            return make_vec3(
                result.x/result.w, result.y/result.w, result.z/result.w
            );
        }
        return triangle_point_[t];
    }
 
    global_index_t triangle_v_global_index(ushort t, index_t llv) const {
        Triangle T = get_triangle(t);
        ushort lv = ushort((llv==0)*T.i + (llv==1)*T.j + (llv==2)*T.k);
        return v_global_index(lv);
    }
 
    index_t triangle_v_local_index(ushort t, index_t llv) const {
        Triangle T = get_triangle(t);
        return index_t((llv==0)*T.i + (llv==1)*T.j + (llv==2)*T.k);
    }
 
    bool triangle_is_user_marked(ushort t) {
        vbw_assert(has_tflags_);
        vbw_assert(t < max_t_);
        return (tflags_[t] != 0);
    }
 
    void triangle_user_mark(ushort t) {
        vbw_assert(has_tflags_);
        vbw_assert(t < max_t_);
        tflags_[t] = 1;
    }
 
    void triangle_user_unmark(ushort t) {
        vbw_assert(has_tflags_);
        vbw_assert(t < max_t_);
        tflags_[t] = 0;
    }
 
    bool cell_has_conflict(const vec4& P) {
        for(
            ushort t = first_triangle();
            t!=END_OF_LIST; t=next_triangle(t)
        ) {
            TriangleWithFlags T = get_triangle_and_flags(t);
            if(triangle_is_in_conflict(T,P)) {
                return true;
            }
        }
        return false;
    }
 
    bool cell_is_totally_in_conflict(const vec4& P) {
        for(
            ushort t = first_triangle();
            t!=END_OF_LIST; t=next_triangle(t)
        ) {
            TriangleWithFlags T = get_triangle_and_flags(t);
            if(!triangle_is_in_conflict(T,P)) {
                return false;
            }
        }
        return true;
    }
 
    index_t triangle_adjacent(index_t t, index_t le) const {
        vbw_assert(t < max_t());
        vbw_assert(le < 3);
        return t_adj_[t][le];
    }
 
 
    void set_triangle_adjacent(index_t t1, index_t le, index_t t2) {
        vbw_assert(t1 < max_t());
        vbw_assert(le < 3);
        vbw_assert(t2 < max_t());
        t_adj_[t1][le] = VBW::ushort(t2);
    }
 
 
 
    index_t triangle_vertex(index_t t, index_t lv) const {
        vbw_assert(t < max_t());
        vbw_assert(lv < 3);
        return t_[t][lv];
    }
 
 
    index_t triangle_find_vertex(index_t t, index_t v) const {
        vbw_assert(t < max_t());
        Triangle T = get_triangle(t);
        index_t result = index_t((T.j == v) + 2*(T.k == v));
        vbw_assert(triangle_vertex(t,result) == v);
        return result;
    }
 
    index_t triangle_find_adjacent(index_t t1, index_t t2) const {
        vbw_assert(t1 < max_t());
        vbw_assert(t2 < max_t());
        Triangle T = t_adj_[t1];
        index_t result = index_t((T.j == t2) + 2*(T.k == t2));
        vbw_assert(triangle_adjacent(t1,result) == t2);
        return result;
    }
 
    bool triangle_is_infinite(index_t t) const {
        vbw_assert(t < max_t());
        Triangle T = get_triangle(t);
        return (
            T.i == VERTEX_AT_INFINITY ||
            T.j == VERTEX_AT_INFINITY ||
            T.k == VERTEX_AT_INFINITY
        );
    }
 
    vec4 vertex_plane(index_t v) const {
        vbw_assert(v < max_v());
        return plane_eqn_[v];
    }
 
    vec3 vertex_plane_normal(index_t v) const {
        vbw_assert(v != VERTEX_AT_INFINITY);
        vbw_assert(v < max_v());
        return make_vec3(
            plane_eqn_[v].x,
            plane_eqn_[v].y,
            plane_eqn_[v].z
        );
    }
 
    bool triangle_is_marked_as_conflict(index_t t) const {
        vbw_assert(t < max_t());
        return (get_triangle_flags(t) & ushort(CONFLICT_MASK)) != 0;
    }
 
    bool triangle_is_in_conflict(
        TriangleWithFlags T, const vec4& eqn
    ) const;
 
    index_t new_triangle(index_t i, index_t j, index_t k) {
        index_t result = first_free_;
        if(result == END_OF_LIST) {
            result = nb_t_;
            ++nb_t_;
            if(nb_t_ > max_t()) {
                grow_t();
            }
        } else {
            first_free_ = index_t(
                get_triangle_flags(first_free_) & ~ushort(CONFLICT_MASK)
            );
        }
        vbw_assert(result < max_t());
        t_[result] = make_triangle_with_flags(
            ushort(i), ushort(j), ushort(k), ushort(first_valid_)
        );
        first_valid_ = result;
        if(has_tflags_) {
            tflags_[result] = 0;
        }
        return result;
    }
 
    index_t new_triangle(
        index_t i, index_t j, index_t k,
        index_t adj0, index_t adj1, index_t adj2
    ) {
        index_t result = new_triangle(i, j, k);
        t_adj_[result] = make_triangle(
            ushort(adj0), ushort(adj1), ushort(adj2)
        );
        return result;
    }
 
    vec4 compute_triangle_point(index_t t) const;
 
    Triangle get_triangle(index_t t) const {
        vbw_assert(t < max_t());
        return t_[t];
    }
 
    ushort get_triangle_flags(index_t t) const {
        vbw_assert(t < max_t());
        return t_[t].flags;
    }
 
    void set_triangle_flags(index_t t, ushort flags) {
        vbw_assert(t < max_t());
        t_[t].flags = flags;
    }
 
    TriangleWithFlags get_triangle_and_flags(index_t t) const {
        vbw_assert(t < max_t());
        return t_[t];
    }
 
    bool triangle_is_marked_as_conflict(index_t t) {
        vbw_assert(t < max_t());
        ushort flg = get_triangle_flags(t);
        return ((flg & ushort(CONFLICT_MASK)) != 0);
    }
 
    index_t max_t() const {
        return max_t_;
    }
 
    index_t max_v() const {
        return max_v_;
    }
 
    void grow_t();
 
    void grow_v();
 
 
    void swap(ConvexCell& other) {
        std::swap(max_t_,other.max_t_);
        std::swap(max_v_,other.max_v_);
        std::swap(t_,other.t_);
        std::swap(t_adj_,other.t_adj_);
        std::swap(plane_eqn_,other.plane_eqn_);
        std::swap(nb_t_,other.nb_t_);
        std::swap(nb_v_,other.nb_v_);
        std::swap(first_free_,other.first_free_);
        std::swap(first_valid_,other.first_valid_);
        std::swap(geometry_dirty_,other.geometry_dirty_);
        std::swap(triangle_point_,other.triangle_point_);
        std::swap(v2t_,other.v2t_);
        std::swap(v2e_,other.v2e_);
        std::swap(vglobal_,other.vglobal_);
        std::swap(has_vglobal_,other.has_vglobal_);
        std::swap(tflags_,other.tflags_);
        std::swap(has_tflags_,other.has_tflags_);
#ifndef STANDALONE_CONVEX_CELL
        std::swap(use_exact_predicates_,other.use_exact_predicates_);
#endif
    }
 
    vec3& stored_triangle_point(ushort t) {
        return triangle_point_[t];
    }
 
    protected:
 
    void connect_triangles();
 
 
    void triangulate_conflict_zone(
        index_t lv, index_t conflict_head, index_t conflict_tail
    );
 
    void set_vertex_plane(index_t v, vec4 P) {
        vbw_assert(v < max_v());
        plane_eqn_[v] = P;
        geometry_dirty_ = true;
    }
 
 
    private:
 
    index_t max_t_;
 
    index_t max_v_;
 
    vector<TriangleWithFlags> t_;
 
    vector<Triangle> t_adj_;
 
    vector<vec4> plane_eqn_;
 
    index_t nb_t_;
 
    index_t nb_v_;
 
    index_t first_free_;
 
    index_t first_valid_;
 
    bool geometry_dirty_;
 
    vector<vec3> triangle_point_;
 
    vector<ushort> v2t_;
 
    vector<uchar> v2e_;
 
    vector<global_index_t> vglobal_;
 
    bool has_vglobal_;
 
    vector<uchar> tflags_;
 
    bool has_tflags_;
 
#ifndef STANDALONE_CONVEX_CELL
    bool use_exact_predicates_;
#endif
 
    friend class GEO::PeriodicDelaunay3d;
    };
}
 
namespace GEO {
    using VBW::ConvexCell;
}
 
#endif