Triangulation.h

Go to the documentation of this file.

/* This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at https://mozilla.org/MPL/2.0/. */
 
#ifndef CDT_vW1vZ0lO8rS4gY4uI4fB
#define CDT_vW1vZ0lO8rS4gY4uI4fB
 
#include "CDTUtils.h"
#include "LocatorKDTree.h"
 
#include <algorithm>
#include <cstdlib>
#include <iterator>
#include <stack>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
namespace CDT
{
 
 
struct CDT_EXPORT VertexInsertionOrder
{
    enum Enum
    {
        Auto,
        AsProvided,
    };
};
 
struct CDT_EXPORT SuperGeometryType
{
    enum Enum
    {
        SuperTriangle, 
        Custom,        
    };
};
 
struct CDT_EXPORT IntersectingConstraintEdges
{
    enum Enum
    {
        NotAllowed, 
        TryResolve, 
        DontCheck,
    };
};
 
typedef unsigned short LayerDepth;
typedef LayerDepth BoundaryOverlapCount;
 
class SourceLocation
{
public:
    SourceLocation(const std::string& file, const std::string& func, int line)
        : m_file(file)
        , m_func(func)
        , m_line(line)
    {}
    const std::string& file() const
    {
        return m_file;
    }
    const std::string& func() const
    {
        return m_func;
    }
    int line() const
    {
        return m_line;
    }
 
private:
    std::string m_file;
    std::string m_func;
    int m_line;
};
 
#define CDT_SOURCE_LOCATION                                                    \
    SourceLocation(std::string(__FILE__), std::string(__func__), __LINE__)
 
class Error : public std::runtime_error
{
public:
    Error(const std::string& description, const SourceLocation& srcLoc)
        : std::runtime_error(
              description + "\nin '" + srcLoc.func() + "' at " + srcLoc.file() +
              ":" + std::to_string(srcLoc.line()))
        , m_description(description)
        , m_srcLoc(srcLoc)
    {}
    const std::string& description() const
    {
        return m_description;
    }
    const SourceLocation& sourceLocation() const
    {
        return m_srcLoc;
    }
 
private:
    std::string m_description;
    SourceLocation m_srcLoc;
};
 
class FinalizedError : public Error
{
public:
    FinalizedError(const SourceLocation& srcLoc)
        : Error(
              "Triangulation was finalized with 'erase...' method. Further "
              "modification is not possible.",
              srcLoc)
    {}
};
 
class DuplicateVertexError : public Error
{
public:
    DuplicateVertexError(
        const VertInd v1,
        const VertInd v2,
        const SourceLocation& srcLoc)
        : Error(
              "Duplicate vertex detected: #" + std::to_string(v1) +
                  " is a duplicate of #" + std::to_string(v2),
              srcLoc)
        , m_v1(v1)
        , m_v2(v2)
    {}
    VertInd v1() const
    {
        return m_v1;
    }
    VertInd v2() const
    {
        return m_v2;
    }
 
private:
    VertInd m_v1, m_v2;
};
 
class IntersectingConstraintsError : public Error
{
public:
    IntersectingConstraintsError(
        const Edge& e1,
        const Edge& e2,
        const SourceLocation& srcLoc)
        : Error(
              "Intersecting constraint edges detected: (" +
                  std::to_string(e1.v1()) + ", " + std::to_string(e1.v2()) +
                  ") intersects (" + std::to_string(e2.v1()) + ", " +
                  std::to_string(e2.v2()) + ")",
              srcLoc)
        , m_e1(e1)
        , m_e2(e2)
    {}
    const Edge& e1() const
    {
        return m_e1;
    }
    const Edge& e2() const
    {
        return m_e2;
    }
 
private:
    Edge m_e1, m_e2;
};
 
 
template <typename T, typename TNearPointLocator = LocatorKDTree<T> >
class CDT_EXPORT Triangulation
{
public:
    typedef std::vector<V2d<T> > V2dVec; 
    V2dVec vertices;                     
    TriangleVec triangles;               
    EdgeUSet fixedEdges; 
 
    unordered_map<Edge, BoundaryOverlapCount> overlapCount;
 
    unordered_map<Edge, EdgeVec> pieceToOriginals;
 
    /*____ API _____*/
    Triangulation();
    explicit Triangulation(VertexInsertionOrder::Enum vertexInsertionOrder);
    Triangulation(
        VertexInsertionOrder::Enum vertexInsertionOrder,
        IntersectingConstraintEdges::Enum intersectingEdgesStrategy,
        T minDistToConstraintEdge);
    Triangulation(
        VertexInsertionOrder::Enum vertexInsertionOrder,
        const TNearPointLocator& nearPtLocator,
        IntersectingConstraintEdges::Enum intersectingEdgesStrategy,
        T minDistToConstraintEdge);
    template <
        typename TVertexIter,
        typename TGetVertexCoordX,
        typename TGetVertexCoordY>
    void insertVertices(
        TVertexIter first,
        TVertexIter last,
        TGetVertexCoordX getX,
        TGetVertexCoordY getY);
    void insertVertices(const std::vector<V2d<T> >& vertices);
    template <
        typename TEdgeIter,
        typename TGetEdgeVertexStart,
        typename TGetEdgeVertexEnd>
    void insertEdges(
        TEdgeIter first,
        TEdgeIter last,
        TGetEdgeVertexStart getStart,
        TGetEdgeVertexEnd getEnd);
    void insertEdges(const std::vector<Edge>& edges);
    template <
        typename TEdgeIter,
        typename TGetEdgeVertexStart,
        typename TGetEdgeVertexEnd>
    void conformToEdges(
        TEdgeIter first,
        TEdgeIter last,
        TGetEdgeVertexStart getStart,
        TGetEdgeVertexEnd getEnd);
    void conformToEdges(const std::vector<Edge>& edges);
    void eraseSuperTriangle();
    void eraseOuterTriangles();
    void eraseOuterTrianglesAndHoles();
    void initializedWithCustomSuperGeometry();
 
    bool isFinalized() const;
 
    std::vector<LayerDepth> calculateTriangleDepths() const;
 
 
    void flipEdge(TriInd iT, TriInd iTopo);
 
    void flipEdge(
        TriInd iT,
        TriInd iTopo,
        VertInd v1,
        VertInd v2,
        VertInd v3,
        VertInd v4,
        TriInd n1,
        TriInd n2,
        TriInd n3,
        TriInd n4);
 
    void removeTriangles(const TriIndUSet& removedTriangles);
 
    TriIndVec& VertTrisInternal();
    const TriIndVec& VertTrisInternal() const;
 
private:
    /*____ Detail __*/
    void addSuperTriangle(const Box2d<T>& box);
    void addNewVertex(const V2d<T>& pos, TriInd iT);
    void insertVertex(VertInd iVert);
    void insertVertex(VertInd iVert, VertInd walkStart);
    void ensureDelaunayByEdgeFlips(VertInd iV1, std::stack<TriInd>& triStack);
    std::vector<Edge> insertVertex_FlipFixedEdges(VertInd iV1);
 
    typedef tuple<IndexSizeType, IndexSizeType, TriInd, TriInd, Index>
        TriangulatePseudoPolygonTask;
 
    void insertEdge(
        Edge edge,
        Edge originalEdge,
        EdgeVec& remaining,
        std::vector<TriangulatePseudoPolygonTask>& tppIterations);
 
    void insertEdgeIteration(
        Edge edge,
        Edge originalEdge,
        EdgeVec& remaining,
        std::vector<TriangulatePseudoPolygonTask>& tppIterations);
 
    typedef tuple<Edge, EdgeVec, BoundaryOverlapCount> ConformToEdgeTask;
 
    void conformToEdge(
        Edge edge,
        EdgeVec originals,
        BoundaryOverlapCount overlaps,
        std::vector<ConformToEdgeTask>& remaining);
 
    void conformToEdgeIteration(
        Edge edge,
        const EdgeVec& originals,
        BoundaryOverlapCount overlaps,
        std::vector<ConformToEdgeTask>& remaining);
 
    tuple<TriInd, VertInd, VertInd> intersectedTriangle(
        VertInd iA,
        const V2d<T>& a,
        const V2d<T>& b,
        T orientationTolerance = T(0)) const;
    std::stack<TriInd> insertVertexInsideTriangle(VertInd v, TriInd iT);
    std::stack<TriInd> insertVertexOnEdge(VertInd v, TriInd iT1, TriInd iT2);
    array<TriInd, 2> trianglesAt(const V2d<T>& pos) const;
    array<TriInd, 2>
    walkingSearchTrianglesAt(VertInd iV, VertInd startVertex) const;
    TriInd walkTriangles(VertInd startVertex, const V2d<T>& pos) const;
    void edgeFlipInfo(
        TriInd iT,
        VertInd iV1,
        TriInd& iTopo,
        VertInd& iV2,
        VertInd& iV3,
        VertInd& iV4,
        TriInd& n1,
        TriInd& n2,
        TriInd& n3,
        TriInd& n4);
    bool isFlipNeeded(VertInd iV1, VertInd iV2, VertInd iV3, VertInd iV4) const;
    void changeNeighbor(TriInd iT, TriInd oldNeighbor, TriInd newNeighbor);
    void changeNeighbor(
        TriInd iT,
        VertInd iVedge1,
        VertInd iVedge2,
        TriInd newNeighbor);
    void triangulatePseudoPolygon(
        const std::vector<VertInd>& poly,
        unordered_map<Edge, TriInd>& outerTris,
        TriInd iT,
        TriInd iN,
        std::vector<TriangulatePseudoPolygonTask>& iterations);
    void triangulatePseudoPolygonIteration(
        const std::vector<VertInd>& poly,
        unordered_map<Edge, TriInd>& outerTris,
        std::vector<TriangulatePseudoPolygonTask>& iterations);
    IndexSizeType findDelaunayPoint(
        const std::vector<VertInd>& poly,
        IndexSizeType iA,
        IndexSizeType iB) const;
    TriInd addTriangle(const Triangle& t); // note: invalidates iterators!
    TriInd addTriangle(); // note: invalidates triangle iterators!
    void finalizeTriangulation(const TriIndUSet& removedTriangles);
    TriIndUSet growToBoundary(std::stack<TriInd> seeds) const;
    void fixEdge(const Edge& edge);
    void fixEdge(const Edge& edge, const Edge& originalEdge);
    void splitFixedEdge(const Edge& edge, const VertInd iSplitVert);
    VertInd addSplitEdgeVertex(
        const V2d<T>& splitVert,
        const TriInd iT,
        const TriInd iTopo);
    VertInd splitFixedEdgeAt(
        const Edge& edge,
        const V2d<T>& splitVert,
        const TriInd iT,
        const TriInd iTopo);
    void makeDummy(TriInd iT);
    void eraseDummies();
    unordered_map<TriInd, LayerDepth> peelLayer(
        std::stack<TriInd> seeds,
        LayerDepth layerDepth,
        std::vector<LayerDepth>& triDepths) const;
 
    void insertVertices_AsProvided(VertInd superGeomVertCount);
    void insertVertices_Randomized(VertInd superGeomVertCount);
    void insertVertices_KDTreeBFS(
        VertInd superGeomVertCount,
        V2d<T> boxMin,
        V2d<T> boxMax);
    std::pair<TriInd, TriInd> edgeTriangles(VertInd a, VertInd b) const;
    bool hasEdge(VertInd a, VertInd b) const;
    void setAdjacentTriangle(const VertInd v, const TriInd t);
    void pivotVertexTriangleCW(VertInd v);
    void tryAddVertexToLocator(const VertInd v);
    void tryInitNearestPointLocator();
 
    std::vector<TriInd> m_dummyTris;
    TNearPointLocator m_nearPtLocator;
    IndexSizeType m_nTargetVerts;
    SuperGeometryType::Enum m_superGeomType;
    VertexInsertionOrder::Enum m_vertexInsertionOrder;
    IntersectingConstraintEdges::Enum m_intersectingEdgesStrategy;
    T m_minDistToConstraintEdge;
    TriIndVec m_vertTris; 
};
 
 
namespace detail
{
 
struct SplitMix64RandGen
{
    typedef unsigned long long uint64;
    uint64 m_state;
    explicit SplitMix64RandGen(uint64 state)
        : m_state(state)
    {}
    explicit SplitMix64RandGen()
        : m_state(0)
    {}
    uint64 operator()()
    {
        uint64 z = (m_state += 0x9e3779b97f4a7c15);
        z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
        z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
        return z ^ (z >> 31);
    }
};
 
template <class RandomIt>
void random_shuffle(RandomIt first, RandomIt last)
{
    detail::SplitMix64RandGen prng;
    typename std::iterator_traits<RandomIt>::difference_type i, n;
    n = last - first;
    for(i = n - 1; i > 0; --i)
    {
        std::swap(first[i], first[prng() % (i + 1)]);
    }
}
 
// backport from c++11
template <class ForwardIt, class T>
void iota(ForwardIt first, ForwardIt last, T value)
{
    while(first != last)
    {
        *first++ = value;
        ++value;
    }
}
 
} // namespace detail
 
//-----------------------
// Triangulation methods
//-----------------------
template <typename T, typename TNearPointLocator>
template <
    typename TVertexIter,
    typename TGetVertexCoordX,
    typename TGetVertexCoordY>
void Triangulation<T, TNearPointLocator>::insertVertices(
    const TVertexIter first,
    const TVertexIter last,
    TGetVertexCoordX getX,
    TGetVertexCoordY getY)
{
    if(isFinalized())
        throw FinalizedError(CDT_SOURCE_LOCATION);
 
    const bool isFirstTime = vertices.empty();
    const T max = std::numeric_limits<T>::max();
    Box2d<T> box = {{max, max}, {-max, -max}};
    if(vertices.empty()) // called first time
    {
        box = envelopBox<T>(first, last, getX, getY);
        addSuperTriangle(box);
    }
    tryInitNearestPointLocator();
 
    const VertInd nExistingVerts = static_cast<VertInd>(vertices.size());
    const VertInd nVerts =
        static_cast<VertInd>(nExistingVerts + std::distance(first, last));
    // optimization, try to pre-allocate tris
    triangles.reserve(triangles.size() + 2 * nVerts);
    vertices.reserve(nVerts);
    m_vertTris.reserve(nVerts);
    for(TVertexIter it = first; it != last; ++it)
        addNewVertex(V2d<T>::make(getX(*it), getY(*it)), noNeighbor);
 
    switch(m_vertexInsertionOrder)
    {
    case VertexInsertionOrder::AsProvided:
        insertVertices_AsProvided(nExistingVerts);
        break;
    case VertexInsertionOrder::Auto:
        isFirstTime ? insertVertices_KDTreeBFS(nExistingVerts, box.min, box.max)
                    : insertVertices_Randomized(nExistingVerts);
        break;
    }
}
 
template <typename T, typename TNearPointLocator>
template <
    typename TEdgeIter,
    typename TGetEdgeVertexStart,
    typename TGetEdgeVertexEnd>
void Triangulation<T, TNearPointLocator>::insertEdges(
    TEdgeIter first,
    const TEdgeIter last,
    TGetEdgeVertexStart getStart,
    TGetEdgeVertexEnd getEnd)
{
    if(isFinalized())
        throw FinalizedError(CDT_SOURCE_LOCATION);
 
    std::vector<TriangulatePseudoPolygonTask> tppIterations;
    EdgeVec remaining;
    for(; first != last; ++first)
    {
        // +3 to account for super-triangle vertices
        const Edge edge(
            VertInd(getStart(*first) + m_nTargetVerts),
            VertInd(getEnd(*first) + m_nTargetVerts));
        insertEdge(edge, edge, remaining, tppIterations);
    }
    eraseDummies();
}
 
template <typename T, typename TNearPointLocator>
template <
    typename TEdgeIter,
    typename TGetEdgeVertexStart,
    typename TGetEdgeVertexEnd>
void Triangulation<T, TNearPointLocator>::conformToEdges(
    TEdgeIter first,
    const TEdgeIter last,
    TGetEdgeVertexStart getStart,
    TGetEdgeVertexEnd getEnd)
{
    if(isFinalized())
        throw FinalizedError(CDT_SOURCE_LOCATION);
 
    tryInitNearestPointLocator();
    // state shared between different runs for performance gains
    std::vector<ConformToEdgeTask> remaining;
    for(; first != last; ++first)
    {
        // +3 to account for super-triangle vertices
        const Edge e(
            VertInd(getStart(*first) + m_nTargetVerts),
            VertInd(getEnd(*first) + m_nTargetVerts));
        conformToEdge(e, EdgeVec(1, e), 0, remaining);
    }
    eraseDummies();
}
 
} // namespace CDT
 
#ifndef CDT_USE_AS_COMPILED_LIBRARY
#include "Triangulation.hpp"
#endif
 
#endif // header-guard