generic_RVD_vertex.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_VERTEX
#define GEOGRAM_VORONOI_GENERIC_RVD_VERTEX
 
#include <geogram/basic/common.h>
#include <geogram/mesh/mesh.h>
#include <geogram/delaunay/delaunay_nn.h>
#include <geogram/basic/assert.h>
#include <geogram/basic/process.h>
#include <geogram/basic/attributes.h>
 
namespace GEOGen {
 
    using GEO::Delaunay;      
    using GEO::index_t;   
    using GEO::signed_index_t;  
    using GEO::coord_index_t;   
    using GEO::Sign;  
    using GEO::Mesh;
 
    template <class T, index_t DIM>
    class small_set {
 
        typedef small_set<T, DIM> thisclass;
 
    public:
        typedef T* iterator;
 
        typedef const T* const_iterator;
 
        typedef T& reference;
 
        typedef T value_type;
 
        small_set() :
            size_(0) {
        }
 
        index_t size() const {
            return size_;
        }
 
        index_t capacity() const {
            return (index_t) DIM;
        }
 
        iterator begin() {
            return data_;
        }
 
        iterator end() {
            return data_ + size_;
        }
 
        iterator end_of_storage() {
            return data_ + DIM;
        }
 
        const_iterator begin() const {
            return data_;
        }
 
        const_iterator end() const {
            return data_ + size_;
        }
 
        const_iterator end_of_storage() const {
            return data_ + DIM;
        }
 
        iterator insert(const T& x) {
            return insert(x, find_i(x));
        }
 
        iterator insert(const T& x, iterator where) {
            if(where == end()) {
                *where = x;
                grow();
                return where;
            }
            if(*where == x) {
                return where;
            }
            grow();
            if(where == end() - 1) {
                *where = x;
                return where;
            }
            for(iterator i = end() - 1; i != where; i--) {
                geo_debug_assert(i != begin());
                *i = *(i - 1);
            }
            *where = x;
#ifdef GEO_DEBUG
            for(iterator i = begin(); i != end() - 1; ++i) {
                geo_debug_assert(*i < *(i + 1));
            }
#endif
            return where;
        }
 
        void clear() {
            size_ = 0;
        }
 
        iterator find(const T& x) {
            iterator result = find_i(x);
            if(*result != x) {
                result = end();
            }
            return result;
        }
 
        const_iterator find(const T& x) const {
            const_iterator result = find_i(x);
            if(*result != x) {
                result = end();
            }
            return result;
        }
 
        void push_back(const T& x) {
#ifdef GEO_DEBUG
            for(iterator i = begin(); i != end(); ++i) {
                geo_debug_assert(*i < x);
            }
#endif
            *end() = x;
            grow();
        }
 
        void print(std::ostream& out) const {
            out << "[ ";
            for(const_iterator it = begin(); it != end(); ++it) {
                out << *it << " ";
            }
            out << "]";
        }
 
        T& operator[] (signed_index_t i) {
            geo_debug_assert(i >= 0);
            geo_debug_assert(begin() + i < end());
            return begin()[i];
        }
 
        const T& operator[] (signed_index_t i) const {
            geo_debug_assert(i >= 0);
            geo_debug_assert(begin() + i < end());
            return begin()[i];
        }
 
    protected:
        void grow() {
            geo_debug_assert(end() != end_of_storage());
            size_++;
        }
 
        // Note: maybe we should start from end() instead of begin()
        // since negative indices are inserted first.
 
        iterator find_i(const T& x) {
            iterator result = begin();
            while(result != end() && *result < x) {
                result++;
            }
            return result;
        }
 
        const_iterator find_i(const T& x) const {
            const_iterator result = begin();
            while(result != end() && *result < x) {
                result++;
            }
            return result;
        }
 
    protected:
        T data_[DIM];
        index_t size_;
    };
 
    template <class T, index_t DIM>
    inline std::ostream& operator<< (
        std::ostream& out,
        const small_set<T, DIM>& S) {
        S.print(out);
        return out;
    }
 
    template <class T, index_t DIM1, index_t DIM2, index_t DIM3>
    inline void sets_intersect(
        const small_set<T, DIM1>& S1,
        const small_set<T, DIM2>& S2,
        small_set<T, DIM3>& I
    ) {
        I.clear();
        auto i1 = S1.begin();
        auto i2 = S2.begin();
        while(i1 < S1.end() && i2 < S2.end()) {
            if(*i1 < *i2) {
                ++i1;
            }
            else if(*i2 < *i1) {
                ++i2;
            }
            else {
                I.push_back(*i1);
                ++i1;
                ++i2;
            }
        }
    }
 
    class SymbolicVertex : public small_set<GEO::signed_index_t, 3> {
 
        typedef SymbolicVertex thisclass;
 
        typedef small_set<GEO::signed_index_t, 3> baseclass;
 
    public:
        SymbolicVertex() :
            v1_(0),
            v2_(0) {
        }
 
        void add_bisector(index_t i) {
            baseclass::insert(signed_index_t(i) + 1);
        }
 
        void add_boundary_facet(index_t i) {
            baseclass::insert(-signed_index_t(i) - 1);
        }
 
        index_t nb_boundary_facets() const {
            index_t result = 0;
            for(auto it = baseclass::begin();
                it != baseclass::end() && *it < 0; ++it) {
                result++;
            }
            return result;
        }
 
        index_t nb_bisectors() const {
            index_t result = 0;
            for(auto it = baseclass::end() - 1;
                it != baseclass::begin() - 1 && *it > 0; --it) {
                result++;
            }
            return result;
        }
 
        static index_t to_unsigned_int(signed_index_t x) {
            geo_debug_assert(x >= 0);
            return (index_t) (x);
        }
 
        index_t bisector(signed_index_t i) const {
            geo_debug_assert(i < (signed_index_t) nb_bisectors());
            return to_unsigned_int((baseclass::end()[-1 - i]) - 1);
        }
 
        index_t boundary_facet(signed_index_t i) const {
            geo_debug_assert(i < (signed_index_t) nb_boundary_facets());
            return to_unsigned_int(-(baseclass::begin()[i]) - 1);
        }
 
        bool has_bisector(index_t i) const {
            return baseclass::find(signed_index_t(i) + 1) != baseclass::end();
        }
 
        bool has_boundary_facet(index_t i) const {
            return baseclass::find(-signed_index_t(i) - 1) != baseclass::end();
        }
 
        index_t get_boundary_vertex() const {
            geo_debug_assert(nb_boundary_facets() == 3);
            geo_debug_assert(v1_ != 0);
            return v1_ - 1;
        }
 
        void get_boundary_edge(index_t& v1, index_t& v2) const {
            geo_debug_assert(nb_boundary_facets() == 2);
            geo_debug_assert(v1_ != 0);
            geo_debug_assert(v2_ != 0);
            v1 = v1_ - 1;
            v2 = v2_ - 1;
        }
 
        void set_boundary_vertex(index_t v) {
            v1_ = v + 1;
            v2_ = 0;
        }
 
        void set_boundary_edge(index_t v1, index_t v2) {
            v1_ = v1 + 1;
            v2_ = v2 + 1;
        }
 
        void copy_boundary_edge_from(const thisclass& rhs) {
            geo_debug_assert(rhs.nb_boundary_facets() == 2);
            geo_debug_assert(rhs.nb_bisectors() == 1);
            geo_debug_assert(rhs.v1_ > 0);
            geo_debug_assert(rhs.v2_ > 0);
            v1_ = rhs.v1_;
            v2_ = rhs.v2_;
        }
 
        bool intersect_symbolic(
            const thisclass& v1,
            const thisclass& v2,
            index_t E
        ) {
 
            // Compute the symbolic representation as the intersection
            // of three planes.
            sets_intersect(v1, v2, *this);
            // this computes the set of planes that contain
            // the edge [v1,v2]
 
            add_bisector(E);   // the intersection is on E.
 
            // Compute the symbolic representation as intersection between
            //    bisector and boundary edge
            // (it's redundant and less elegant than the representation
            //  as planes interactions,
            //  but we need this to handle degenerate configurations properly,
            //  and to use exact predicates with original boundary vertices
            //  coordinates).
 
            if(nb_boundary_facets() == 2) {
                // If *this is on the intersection of two boundary facets,
                // then *this is on
                // a boundary edge, and we need to retrieve the indices of the
                // two extremities of this boundary edge.
 
                index_t nb1 = v1.nb_boundary_facets();
                index_t nb2 = v2.nb_boundary_facets();
                if(nb1 == 3 && nb2 == 3) {
                    // If v1 and v2 are boundary vertices,
                    // then I is on the boundary
                    // edge that connects v1 and v2
                    set_boundary_edge(
                        v1.get_boundary_vertex(),
                        v2.get_boundary_vertex()
                    );
                } else if(nb1 == 2) {
                    geo_debug_assert(nb_boundary_facets() == 2);
                    // If v1 is on a boundary edge,
                    // then I is on the same boundary edge as v1
                    copy_boundary_edge_from(v1);
                } else if(nb2 == 2) {
                    geo_debug_assert(nb_boundary_facets() == 2);
                    // If v2 is on a boundary edge,
                    // then I is on the same boundary edge as v2
                    copy_boundary_edge_from(v2);
                }
            }
 
            // Sanity check: problem detected here, we
            // notify the caller that will use a workaround
            // (see clip_by_plane())
            if(baseclass::size() != 3) {
                return false;
            }
            return true;
        }
 
    private:
        index_t v1_;
        index_t v2_;
    };
 
 
    class PointAllocator {
    public:
        PointAllocator(coord_index_t dim) :
            size_(0),
            capacity_(0),
            dimension_(dim) {
        }
 
        double* new_item() {
            if(size_ == capacity_) {
                grow();
            }
            size_++;
            return item(size_ - 1);
        }
 
        void clear() {
            size_ = 0;
        }
 
        ~PointAllocator() {
            for(index_t c = 0; c < chunks_.size(); c++) {
                GEO::Memory::aligned_free(chunks_[c]);
            }
        }
 
        coord_index_t dimension() const {
            return dimension_;
        }
 
    protected:
        enum {
            CHUNK_SHIFT = 8,
            CHUNK_SIZE = 1 << CHUNK_SHIFT,
            CHUNK_MASK = CHUNK_SIZE - 1
        };
 
        void grow() {
            chunks_.push_back(
                reinterpret_cast<double*>(
                    GEO::Memory::aligned_malloc(
                        index_t(CHUNK_SIZE) * dimension_ * sizeof(double)
                    )
                )
            );
            capacity_ += CHUNK_SIZE;
        }
 
        double* item(index_t i) {
            geo_debug_assert(i < size_);
            return &(chunks_[i >> CHUNK_SHIFT][(i & CHUNK_MASK) * dimension_]);
        }
 
    private:
        index_t size_;
        index_t capacity_;
        coord_index_t dimension_;
        std::vector<double*> chunks_;
    };
 
    enum {
        ORIGINAL = 1,  
        INTERSECT = 2  
    };
 
    typedef index_t EdgeFlags;
 
    typedef index_t EdgeFlag;
 
    class Vertex {
 
        typedef Vertex thisclass;
 
    public:
        Vertex(
            const double* p, double w, signed_index_t f,
            const SymbolicVertex& sym
        ) :
            point_(p),
            weight_(w),
            f_(f),
            seed_(-1),
            sym_(sym),
            flags_(ORIGINAL) {
        }
 
        Vertex(const double* p, double w, signed_index_t f) :
            point_(p),
            weight_(w),
            f_(f),
            seed_(-1),
            flags_(ORIGINAL) {
        }
 
        Vertex() :
            point_(nullptr),
            weight_(1.0),
            f_(-1),
            seed_(-1),
            flags_(0) {
        }
 
        const double* point() const {
            return point_;
        }
 
        void set_point(const double* p) {
            point_ = p;
        }
 
        double weight() const {
            return weight_;
        }
 
        void set_weight(double w) {
            weight_ = w;
        }
 
        signed_index_t adjacent_seed() const {
            return seed_;
        }
 
        void set_adjacent_seed(signed_index_t s) {
            seed_ = s;
        }
 
        const SymbolicVertex& sym() const {
            return sym_;
        }
 
        SymbolicVertex& sym() {
            return sym_;
        }
 
        signed_index_t adjacent_facet() const {
            return f_;
        }
 
        void set_adjacent_facet(signed_index_t f) {
            f_ = f;
        }
 
        operator const double* () const {
            return point_;
        }
 
        void clear() {
            flags_ = 0;
            f_ = -1;
        }
 
        void set_flag(EdgeFlag f) {
            flags_ |= f;
        }
 
        void unset_flag(EdgeFlag f) {
            flags_ &= ~f;
        }
 
        bool check_flag(EdgeFlag f) const {
            return (flags_ & f) != 0;
        }
 
        void copy_edge_from(const Vertex& rhs) {
            set_adjacent_facet(rhs.adjacent_facet());
            flags_ = rhs.flags_;
        }
 
        template <index_t DIM>
        void intersect_geom(
            PointAllocator& target_intersections,
            const Vertex& vq1, const Vertex& vq2,
            const double* p1, const double* p2
        ) {
            const double* q1 = vq1.point();
            const double* q2 = vq2.point();
            double* Ipoint = target_intersections.new_item();
            set_point(Ipoint);
            double d = 0.0, l1 = 0.0, l2 = 0.0;
            for(coord_index_t c = 0; c < DIM; ++c) {
                double n = p1[c] - p2[c];
                d -= n * (p2[c] + p1[c]);
                l1 += q2[c] * n;
                l2 += q1[c] * n;
            }
            d = 0.5 * d;
            l1 = ::fabs(l1 + d);
            l2 = ::fabs(l2 + d);
            double l12 = l1 + l2;
            if(l12 > 1e-30) {
                l1 /= l12;
                l2 /= l12;
            } else {
                l1 = 0.5;
                l2 = 0.5;
            }
            for(coord_index_t c = 0; c < DIM; ++c) {
                Ipoint[c] = l1 * q1[c] + l2 * q2[c];
            }
            set_weight(l1 * vq1.weight() + l2 * vq2.weight());
        }
 
        template <index_t DIM>
        Sign side_fast(
            const double* p1, const double* p2
        ) const {
            double r = 0.0;
            for(index_t c = 0; c < DIM; ++c) {
                r += GEO::geo_sqr(p2[c] - point()[c]);
                r -= GEO::geo_sqr(p1[c] - point()[c]);
            }
            return GEO::geo_sgn(r);
        }
 
    private:
        const double* point_;
        double weight_;
 
        signed_index_t f_;
 
        signed_index_t seed_;
 
        SymbolicVertex sym_;
 
        EdgeFlags flags_;
    };
}
 
/*
  namespace GEO {
  template <> struct can_be_used_as_attribute<GEOGen::SymbolicVertex> {
  static constexpr auto value = std::integral_constant<bool,true>();
  };
  }
*/
 
 
#endif