CDT_2d.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_DELAUNAY_CDT_2D
#define GEOGRAM_DELAUNAY_CDT_2D
 
#include <geogram/basic/common.h>
#include <geogram/basic/geometry.h>
#include <geogram/numerics/predicates.h>
#include <geogram/numerics/exact_geometry.h>
#include <geogram/mesh/index.h>
#include <functional>
 
namespace GEO {
 
    struct CDT2d_ConstraintWalker;
 
    class GEOGRAM_API CDTBase2d {
    public:
    CDTBase2d();
 
    virtual ~CDTBase2d();
 
    virtual void clear();
 
    void insert_constraint(index_t i, index_t j);
 
    void remove_external_triangles(
        bool remove_internal_holes=false
    );
 
    void set_delaunay(bool delaunay) {
        delaunay_ = delaunay;
    }
 
    index_t nT() const {
        return T_.size()/3;
    }
 
    index_t nv() const {
        return nv_;
    }
 
    index_t ncnstr() const {
        return ncnstr_;
    }
 
    index_t Tv(index_t t, index_t lv) const {
        geo_debug_assert(t<nT());
        geo_debug_assert(lv<3);
        return T_[3*t+lv];
    }
 
    index_t Tv_find(index_t t, index_t v) const {
        geo_debug_assert(t<nT());
        geo_debug_assert(v<nv());
        return find_3(T_.data()+3*t, v);
    }
 
    index_t Tadj(index_t t, index_t le) const {
        geo_debug_assert(t<nT());
        geo_debug_assert(le<3);
        return Tadj_[3*t+le];
    }
 
    index_t Tadj_find(index_t t1, index_t t2) const {
        geo_debug_assert(t1<nT());
        geo_debug_assert(t2<nT());
        return find_3(Tadj_.data()+3*t1, t2);
    }
 
    index_t vT(index_t v) const {
        geo_debug_assert(v < nv());
        return v2T_[v];
    }
 
 
    index_t Tedge_cnstr_first(index_t t, index_t le) const {
        geo_debug_assert(t < nT());
        geo_debug_assert(le < 3);
        return Tecnstr_first_[3*t+le];
    }
 
    index_t edge_cnstr_next(index_t ecit) const {
        return ecnstr_next_[ecit];
    }
 
    index_t edge_cnstr(index_t ecit) const {
        return ecnstr_val_[ecit];
    }
 
    index_t Tedge_cnstr_nb(index_t t, index_t le) const {
        index_t result = 0;
        for(
            index_t ecit = Tedge_cnstr_first(t,le);
            ecit != index_t(-1);
            ecit = edge_cnstr_next(ecit)
        ) {
            ++result;
        }
        return result;
    }
 
 
    virtual void save(const std::string& filename) const = 0;
 
    bool Tedge_is_Delaunay(index_t t, index_t le) const;
 
    protected:
    virtual void begin_insert_transaction();
    virtual void commit_insert_transaction();
    virtual void rollback_insert_transaction();
 
    index_t insert(index_t v, index_t hint = index_t(-1));
 
    void create_enclosing_triangle(index_t v1, index_t v2, index_t v3);
 
    void create_enclosing_quad(
        index_t v1, index_t v2, index_t v3, index_t v4
    );
 
    void Tset_flag(index_t t, index_t flag) {
        geo_debug_assert(t < nT());
        geo_debug_assert(flag < 8);
        Tflags_[t] |= Numeric::uint8(1u << flag);
    }
 
    void Treset_flag(index_t t, index_t flag) {
        geo_debug_assert(t < nT());
        geo_debug_assert(flag < 8);
        Tflags_[t] &= Numeric::uint8(~(1u << flag));
    }
 
    bool Tflag_is_set(index_t t, index_t flag) {
        geo_debug_assert(t < nT());
        geo_debug_assert(flag < 8);
        return ((Tflags_[t] & (1u << flag)) != 0);
    }
 
    enum {
        DLIST_S_ID=0,
        DLIST_Q_ID=1,
        DLIST_N_ID=2,
        DLIST_NB=3
    };
 
    enum {
        T_MARKED_FLAG  = DLIST_NB,
        T_VISITED_FLAG = DLIST_NB+1
    };
 
    bool Tis_in_list(index_t t) const {
        return (
            (Tflags_[t] &
             Numeric::uint8((1 << DLIST_NB)-1)
            ) != 0
        );
    }
 
    struct DList {
        DList(CDTBase2d& cdt, index_t list_id) :
            cdt_(cdt), list_id_(list_id),
            back_(index_t(-1)), front_(index_t(-1)) {
            geo_debug_assert(list_id < DLIST_NB);
        }
 
        DList(CDTBase2d& cdt) :
            cdt_(cdt), list_id_(index_t(-1)),
            back_(index_t(-1)), front_(index_t(-1)) {
        }
 
        void initialize(index_t list_id) {
            geo_debug_assert(!initialized());
            geo_debug_assert(list_id < DLIST_NB);
            list_id_ = list_id;
        }
 
        bool initialized() const {
            return (list_id_ != index_t(-1));
        }
 
        ~DList() {
            if(initialized()) {
                clear();
            }
        }
 
        bool empty() const {
            geo_debug_assert(initialized());
            geo_debug_assert(
                (back_==index_t(-1))==(front_==index_t(-1))
            );
            return (back_==index_t(-1));
        }
 
        bool contains(index_t t) const {
            geo_debug_assert(initialized());
            return cdt_.Tflag_is_set(t, list_id_);
        }
 
        index_t front() const {
            geo_debug_assert(initialized());
            return front_;
        }
 
        index_t back() const {
            geo_debug_assert(initialized());
            return back_;
        }
 
        index_t next(index_t t) const {
            geo_debug_assert(initialized());
            geo_debug_assert(contains(t));
            return cdt_.Tnext_[t];
        }
 
        index_t prev(index_t t) const {
            geo_debug_assert(initialized());
            geo_debug_assert(contains(t));
            return cdt_.Tprev_[t];
        }
 
        void clear() {
            for(index_t t=front_; t!=index_t(-1); t = cdt_.Tnext_[t]) {
                cdt_.Treset_flag(t,list_id_);
            }
            back_ = index_t(-1);
            front_ = index_t(-1);
        }
 
        index_t size() const {
            geo_debug_assert(initialized());
            index_t result = 0;
            for(index_t t=front(); t!=index_t(-1); t = next(t)) {
                ++result;
            }
            return result;
        }
 
        void push_back(index_t t) {
            geo_debug_assert(initialized());
            geo_debug_assert(!cdt_.Tis_in_list(t));
            cdt_.Tset_flag(t,list_id_);
            if(empty()) {
                back_ = t;
                front_ = t;
                cdt_.Tnext_[t] = index_t(-1);
                cdt_.Tprev_[t] = index_t(-1);
            } else {
                cdt_.Tnext_[t] = index_t(-1);
                cdt_.Tnext_[back_] = t;
                cdt_.Tprev_[t] = back_;
                back_ = t;
            }
        }
 
        index_t pop_back() {
            geo_debug_assert(initialized());
            geo_debug_assert(!empty());
            index_t t = back_;
            back_ = cdt_.Tprev_[back_];
            if(back_ == index_t(-1)) {
                geo_debug_assert(front_ == t);
                front_ = index_t(-1);
            } else {
                cdt_.Tnext_[back_] = index_t(-1);
            }
            geo_debug_assert(contains(t));
            cdt_.Treset_flag(t,list_id_);
            return t;
        }
 
        void push_front(index_t t) {
            geo_debug_assert(initialized());
            geo_debug_assert(!cdt_.Tis_in_list(t));
            cdt_.Tset_flag(t,list_id_);
            if(empty()) {
                back_ = t;
                front_ = t;
                cdt_.Tnext_[t] = index_t(-1);
                cdt_.Tprev_[t] = index_t(-1);
            } else {
                cdt_.Tprev_[t] = index_t(-1);
                cdt_.Tprev_[front_] = t;
                cdt_.Tnext_[t] = front_;
                front_ = t;
            }
        }
 
        index_t pop_front() {
            geo_debug_assert(initialized());
            geo_debug_assert(!empty());
            index_t t = front_;
            front_ = cdt_.Tnext_[front_];
            if(front_ == index_t(-1)) {
                geo_debug_assert(back_ == t);
                back_ = index_t(-1);
            } else {
                cdt_.Tprev_[front_] = index_t(-1);
            }
            geo_debug_assert(contains(t));
            cdt_.Treset_flag(t,list_id_);
            return t;
        }
 
        void remove(index_t t) {
            geo_debug_assert(initialized());
            if(t == front_) {
                pop_front();
            } else if(t == back_) {
                pop_back();
            } else {
                geo_debug_assert(contains(t));
                index_t t_prev = cdt_.Tprev_[t];
                index_t t_next = cdt_.Tnext_[t];
                cdt_.Tprev_[t_next] = t_prev;
                cdt_.Tnext_[t_prev] = t_next;
                cdt_.Treset_flag(t,list_id_);
            }
        }
 
        void show(std::ostream& out = std::cerr) const {
            switch(list_id_) {
            case DLIST_S_ID:
                out << "S";
                break;
            case DLIST_Q_ID:
                out << "Q";
                break;
            case DLIST_N_ID:
                out << "N";
                break;
            case index_t(-1):
                out << "<uninitialized list>";
                break;
            default:
                out << "<unknown list id:" << list_id_ << ">";
                break;
            }
            out << "=";
            for(index_t t=front(); t!=index_t(-1); t = next(t)) {
                out << t << ";";
            }
            out << std::endl;
        }
 
    private:
        CDTBase2d& cdt_;
        index_t list_id_;
        index_t back_;
        index_t front_;
    };
 
    void insert_vertex_in_edge(index_t v, index_t t, index_t le, DList& S);
 
    void insert_vertex_in_edge(index_t v, index_t t, index_t le) {
        DList S(*this);
        insert_vertex_in_edge(v,t,le,S);
    }
 
    void insert_vertex_in_triangle(index_t v, index_t t, DList& S);
 
    index_t find_intersected_edges(index_t i, index_t j, DList& Q);
 
    void walk_constraint_v(CDT2d_ConstraintWalker& W);
 
    void walk_constraint_t(CDT2d_ConstraintWalker& W, DList& Q);
 
    void constrain_edges(index_t i, index_t j, DList& Q, DList& N);
 
    void Delaunayize_vertex_neighbors(index_t from_v);
 
    void Delaunayize_vertex_neighbors(index_t v, DList& S);
 
    void Delaunayize_new_edges(DList& N);
 
 
    void Tset(
        index_t t,
        index_t v1,   index_t v2,   index_t v3,
        index_t adj1, index_t adj2, index_t adj3,
        index_t e1cnstr = index_t(-1),
        index_t e2cnstr = index_t(-1),
        index_t e3cnstr = index_t(-1)
    ) {
        geo_debug_assert(t < nT());
        geo_debug_assert(v1 < nv());
        geo_debug_assert(v2 < nv());
        geo_debug_assert(v3 < nv());
        geo_debug_assert(adj1 < nT() || adj1 == index_t(-1));
        geo_debug_assert(adj2 < nT() || adj2 == index_t(-1));
        geo_debug_assert(adj3 < nT() || adj3 == index_t(-1));
        geo_debug_assert(v1 != v2);
        geo_debug_assert(v2 != v3);
        geo_debug_assert(v3 != v1);
        geo_debug_assert(adj1 != adj2 || adj1 == index_t(-1));
        geo_debug_assert(adj2 != adj3 || adj2 == index_t(-1));
        geo_debug_assert(adj3 != adj1 || adj3 == index_t(-1));
        geo_debug_assert(orient2d(v1,v2,v3) != ZERO);
        T_[3*t  ]    = v1;
        T_[3*t+1]    = v2;
        T_[3*t+2]    = v3;
        Tadj_[3*t  ] = adj1;
        Tadj_[3*t+1] = adj2;
        Tadj_[3*t+2] = adj3;
        Tecnstr_first_[3*t]   = e1cnstr;
        Tecnstr_first_[3*t+1] = e2cnstr;
        Tecnstr_first_[3*t+2] = e3cnstr;
        v2T_[v1] = t;
        v2T_[v2] = t;
        v2T_[v3] = t;
    }
 
    void Trot(index_t t, index_t lv) {
        geo_debug_assert(t < nT());
        geo_debug_assert(lv < 3);
        if(lv != 0) {
            index_t i = 3*t+lv;
            index_t j = 3*t+((lv+1)%3);
            index_t k = 3*t+((lv+2)%3);
            Tset(
                t,
                T_[i], T_[j], T_[k],
                Tadj_[i], Tadj_[j], Tadj_[k],
                Tecnstr_first_[i], Tecnstr_first_[j], Tecnstr_first_[k]
            );
        }
    }
 
    void swap_edge(index_t t1, bool swap_t1_t2=false);
 
    void Tadj_set(index_t t, index_t le, index_t adj) {
        geo_debug_assert(t < nT());
        geo_debug_assert(adj < nT());
        geo_debug_assert(le < 3);
        Tadj_[3*t+le] = adj;
    }
 
    index_t Topp(index_t t, index_t e=0) const {
        index_t t2 = Tadj(t,e);
        if(t2 == index_t(-1)) {
            return index_t(-1);
        }
        index_t e2 = Tadj_find(t2,t);
        return Tv(t2,e2);
    }
 
    void Tadj_back_connect(
        index_t t1, index_t le1, index_t prev_t2_adj_e2
    ) {
        geo_debug_assert(t1 < nT());
        geo_debug_assert(le1 < 3);
        index_t t2 = Tadj(t1,le1);
        if(t2 == index_t(-1)) {
            return;
        }
        index_t le2 = Tadj_find(t2,prev_t2_adj_e2);
        Tadj_set(t2,le2,t1);
        Tset_edge_cnstr_first(t1,le1,Tedge_cnstr_first(t2,le2));
    }
 
    index_t Tnew() {
        index_t t = nT();
        index_t nc = (t+1)*3; // new number of corners
        T_.resize(nc, index_t(-1));
        Tadj_.resize(nc, index_t(-1));
        Tecnstr_first_.resize(nc, index_t(-1));
        Tflags_.resize(t+1,0);
        Tnext_.resize(t+1,index_t(-1));
        Tprev_.resize(t+1,index_t(-1));
        return t;
    }
 
    void Tset_edge_cnstr_first(
        index_t t, index_t le, index_t ecit
    ) {
        geo_debug_assert(t < nT());
        geo_debug_assert(le < 3);
        Tecnstr_first_[3*t+le] = ecit;
    }
 
    void Tadd_edge_cnstr(
        index_t t, index_t le, index_t cnstr_id
    ) {
        geo_debug_assert(t < nT());
        geo_debug_assert(le < 3);
        // Check whether the edge is already constrained with the
        // same constraint.
        // TODO (if possible): understand how this can happen and
        // remove this bloc of code that is not super elegant
        // (it seems to be when we arrive at j and coming from a vertex
        // traversed by the edge, both conditions make the constraint
        // added to the traversed edge).
        for(
            index_t ecit = Tedge_cnstr_first(t,le);
            ecit != index_t(-1);
            ecit = edge_cnstr_next(ecit)
        ) {
            if(edge_cnstr(ecit) == cnstr_id) {
                return;
            }
        }
        ecnstr_val_.push_back(cnstr_id);
        ecnstr_next_.push_back(Tedge_cnstr_first(t,le));
        Tset_edge_cnstr_first(t,le, ecnstr_val_.size()-1);
    }
 
    void Tadd_edge_cnstr_with_neighbor(
        index_t t, index_t le, index_t cnstr_id
    ) {
        geo_debug_assert(t < nT());
        geo_debug_assert(le < 3);
#ifdef GEO_DEBUG
        index_t t_e_cnstr_first = Tedge_cnstr_first(t,le);
#endif
        Tadd_edge_cnstr(t, le, cnstr_id);
        index_t t2 = Tadj(t,le);
        if(t2 != index_t(-1)) {
            index_t le2 = Tadj_find(t2,t);
            // Sanity check: make sure the two edges always share the
            // same constraint list.
            geo_debug_assert(Tedge_cnstr_first(t2,le2) == t_e_cnstr_first);
            Tset_edge_cnstr_first(t2,le2,Tedge_cnstr_first(t,le));
        }
    }
 
    bool Tedge_is_constrained(index_t t, index_t le) const {
        return (Tedge_cnstr_first(t,le) != index_t(-1));
    }
 
    void for_each_T_around_v(
        index_t v, std::function<bool(index_t t, index_t lv)> doit
    ) {
        index_t t = vT(v);
        index_t lv = index_t(-1);
        do {
            lv = Tv_find(t,v);
            if(doit(t,lv)) {
                return;
            }
            t = Tadj(t, (lv+1)%3);
        } while(t != vT(v) && t != index_t(-1));
 
        // We are done, this was an interior vertex
        if(t != index_t(-1)) {
            return;
        }
 
        // It was a vertex on the border, so we need
        // to traverse the triangle fan in the other
        // direction until we reach the border again
        t = vT(v);
        lv = Tv_find(t,v);
        t = Tadj(t, (lv+2)%3);
        while(t != index_t(-1)) {
            lv = Tv_find(t,v);
            if(doit(t,lv)) {
                return;
            }
            t = Tadj(t, (lv+2)%3);
        }
    }
 
 
    index_t locate(
        index_t v, index_t hint = index_t(-1), Sign* orient = nullptr
    ) const;
 
    bool is_convex_quad(index_t t) const;
 
    virtual Sign orient2d(index_t i,index_t j,index_t k) const=0;
 
    virtual Sign incircle(index_t i,index_t j,index_t k,index_t l) const=0;
 
    virtual index_t create_intersection(
        index_t E1, index_t i, index_t j,
        index_t E2, index_t k, index_t l
    ) = 0;
 
    static inline index_t find_3(const index_t* T, index_t v) {
        // The following expression is 10% faster than using
        // if() statements. This uses the C++ norm, that
        // ensures that the 'true' boolean value converted to
        // an int is always 1. With most compilers, this avoids
        // generating branching instructions.
        // Thank to Laurent Alonso for this idea.
        index_t result = index_t( (T[1] == v) | ((T[2] == v) * 2) );
        // Sanity check, important if it was T[0], not explicitly
        // tested (detects input that does not meet the precondition).
        geo_debug_assert(T[result] == v);
        return result;
    }
 
    /*******************************************************************/
 
    void remove_marked_triangles();
 
    /******************** Debugging ************************************/
 
    void Tcheck(index_t t) const {
        if(t == index_t(-1)) {
            return;
        }
        for(index_t e=0; e<3; ++e) {
            geo_assert(Tv(t,e) != Tv(t,(e+1)%3));
            if(Tadj(t,e) == index_t(-1)) {
                continue;
            }
            geo_assert(Tadj(t,e) != Tadj(t,(e+1)%3));
            index_t t2 = Tadj(t,e);
            index_t e2 = Tadj_find(t2,t);
            geo_assert(Tadj(t2,e2) == t);
        }
    }
 
    void debug_Tcheck(index_t t) const {
#ifdef GEO_DEBUG
        Tcheck(t);
#else
        geo_argused(t);
#endif
    }
 
    void check_combinatorics() const {
        for(index_t t=0; t<nT(); ++t) {
            Tcheck(t);
        }
    }
 
    void debug_check_combinatorics() const {
#ifdef GEO_DEBUG
        check_combinatorics();
#endif
    }
 
    virtual void check_geometry() const;
 
    void debug_check_geometry() const {
#ifdef GEO_DEBUG
        check_geometry();
#endif
    }
 
 
    public:
    void check_consistency() const {
        check_combinatorics();
        check_geometry();
    }
 
    protected:
    void debug_check_consistency() const {
        debug_check_combinatorics();
        debug_check_geometry();
    }
 
    bool segment_segment_intersect(
        index_t u1, index_t u2, index_t v1, index_t v2
    ) const {
        if(orient2d(u1,u2,v1)*orient2d(u1,u2,v2) > 0) {
            return false;
        }
        return (orient2d(v1,v2,u1)*orient2d(v1,v2,u2) < 0);
    }
 
    bool segment_edge_intersect(
        index_t v1, index_t v2, index_t t, index_t le
    ) const {
        index_t u1 = Tv(t,(le + 1)%3);
        index_t u2 = Tv(t,(le + 2)%3);
        return segment_segment_intersect(u1,u2,v1,v2);
    }
 
    void check_edge_intersections(
        index_t v1, index_t v2, const DList& Q
    );
 
    typedef std::pair<index_t, index_t> Edge;
 
    index_t eT(Edge E) {
        index_t v1 = E.first;
        index_t v2 = E.second;
        index_t result = index_t(-1);
        for_each_T_around_v(
            v1, [&](index_t t, index_t lv)->bool {
                if(Tv(t, (lv+1)%3) == v2) {
                    if(Tv(t, (lv+2)%3) != v1) {
                        Trot(t, (lv+2)%3);
                    }
                    result = t;
                    return true;
                } else if(Tv(t, (lv+1)%3) == v1) {
                    if(Tv(t, (lv+2)%3) != v2) {
                        Trot(t, (lv+2)%3);
                    }
                    result = t;
                    return true;
                }
                return false;
            }
        );
        geo_debug_assert(result != index_t(-1));
        geo_debug_assert(
            (Tv(result,1) == v1 && Tv(result,2) == v2) ||
            (Tv(result,1) == v2 && Tv(result,2) == v1)
        );
        return result;
    }
 
    index_t locate_naive(
        index_t v, index_t hint = index_t(-1), Sign* orient = nullptr
    ) const;
 
    void constrain_edges_naive(
        index_t i, index_t j, DList& Q, vector<Edge>& N
    );
 
    void Delaunayize_new_edges_naive(vector<Edge>& N);
 
    protected:
    index_t nv_;
    index_t ncnstr_;
    vector<index_t> T_;        
    vector<index_t> Tadj_;     
    vector<index_t> v2T_;      
    vector<uint8_t> Tflags_;   
    vector<index_t> Tecnstr_first_;  
    vector<index_t> ecnstr_val_;     
    vector<index_t> ecnstr_next_;    
    vector<index_t> Tnext_;    
    vector<index_t> Tprev_;    
    bool delaunay_;            
    Sign orient_012_;          
    bool exact_incircle_;      
    bool exact_intersections_; 
    };
 
    /*****************************************************************/
 
    class GEOGRAM_API CDT2d: public CDTBase2d {
    public:
 
        CDT2d();
 
        ~CDT2d() override;
 
        void clear() override;
 
        void create_enclosing_triangle(
            const vec2& p1, const vec2& p2, const vec2& p3
        );
 
        void create_enclosing_quad(
            const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4
        );
 
 
        void create_enclosing_rectangle(
            double x1, double y1, double x2, double y2
        ) {
            create_enclosing_quad(
                vec2(x1,y1),
                vec2(x2,y1),
                vec2(x2,y2),
                vec2(x1,y2)
            );
        }
 
        index_t insert(const vec2& p, index_t hint = index_t(-1)) {
            debug_check_consistency();
            point_.push_back(p);
            index_t v = CDTBase2d::insert(point_.size()-1, hint);
            // If inserted point already existed in
            // triangulation, then nv() did not increase
            if(point_.size() > nv()) {
                point_.pop_back();
            }
            debug_check_consistency();
            return v;
        }
 
        void insert(
            index_t nb_points, const double* points,
            index_t* indices = nullptr,
            bool remove_unreferenced_vertices = false
        );
 
        void save(const std::string& filename) const override;
 
        const vec2 point(index_t v) const {
            geo_debug_assert(v < nv());
            return point_[v];
        }
 
    protected:
        Sign orient2d(index_t i, index_t j, index_t k) const override;
 
        Sign incircle(index_t i,index_t j,index_t k,index_t l) const override;
 
        index_t create_intersection(
            index_t E1, index_t i, index_t j,
            index_t E2, index_t k, index_t l
        ) override;
 
    protected:
        vector<vec2> point_;
    };
 
    /*****************************************************************/
 
    class GEOGRAM_API ExactCDT2d : public CDTBase2d {
    public:
        typedef exact::vec2h ExactPoint;
 
        ExactCDT2d();
 
        ~ExactCDT2d() override;
 
        void clear() override;
 
        index_t insert(
            const ExactPoint& p, index_t id=0, index_t hint = index_t(-1)
        );
 
        void insert_constraint(index_t v1, index_t v2, index_t operand_bits=0) {
            constraints_.push_back(bindex(v1,v2,bindex::KEEP_ORDER));
            cnstr_operand_bits_.push_back(operand_bits);
            CDTBase2d::insert_constraint(v1,v2);
        }
 
        void create_enclosing_quad(
            const ExactPoint& p1, const ExactPoint& p2,
            const ExactPoint& p3, const ExactPoint& p4
        );
 
        void create_enclosing_rectangle(
            double x1, double y1, double x2, double y2
        ) {
            create_enclosing_quad(
                ExactPoint(vec2(x1,y1)),
                ExactPoint(vec2(x2,y1)),
                ExactPoint(vec2(x2,y2)),
                ExactPoint(vec2(x1,y2))
            );
        }
 
        const ExactPoint& point(index_t v) const {
            geo_debug_assert(v < nv());
            return point_[v];
        }
 
        index_t vertex_id(index_t v) const {
            geo_debug_assert(v < nv());
            return id_[v];
        }
 
        void set_vertex_id(index_t v, index_t id) {
            geo_debug_assert(v < nv());
            id_[v] = id;
        }
 
        void classify_triangles(
            const std::string& boolean_expression, bool mark_only=false
        );
 
        void save(const std::string& filename) const override;
 
    protected:
        void add_point(const ExactPoint& p, index_t id = index_t(-1));
        void begin_insert_transaction() override;
        void commit_insert_transaction() override;
        void rollback_insert_transaction() override;
 
        Sign orient2d(index_t i, index_t j, index_t k) const override;
 
        Sign incircle(index_t i,index_t j,index_t k,index_t l) const override;
 
        index_t create_intersection(
            index_t E1, index_t i, index_t j,
            index_t E2, index_t k, index_t l
        ) override;
 
    protected:
        vector<ExactPoint> point_;
#ifndef GEOGRAM_USE_EXACT_NT
        vector<double> length_;
#endif
        vector<index_t> id_;
        vector<index_t> cnstr_operand_bits_;
        vector<index_t> facet_inclusion_bits_;
        mutable std::map<trindex, Sign> pred_cache_;
        bool use_pred_cache_insert_buffer_;
        mutable std::vector<std::pair<trindex, Sign>> pred_cache_insert_buffer_;
        vector<bindex> constraints_;
    };
 
    /*****************************************************************/
 
}
 
#endif