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 CDT_EXPORT Error : public std::runtime_error
{
public:
    Error(const std::string& description, const SourceLocation& srcLoc)
        : std::runtime_error(
              description + "\nin '" + srcLoc.func() + "' at " + srcLoc.file() +
              ":" + CDT::to_string(srcLoc.line()))
        , m_description(description)
        , m_srcLoc(srcLoc)
    {}

    virtual ~Error() CDT_NOEXCEPT
    {}

    const std::string& description() const
    {
        return m_description;
    }

    const SourceLocation& sourceLocation() const
    {
        return m_srcLoc;
    }
 
private:
    std::string m_description;
    SourceLocation m_srcLoc;
};

class CDT_EXPORT FinalizedError : public Error
{
public:
    FinalizedError(const SourceLocation& srcLoc)
        : Error(
              "Triangulation was finalized with 'erase...' method. Further "
              "modification is not possible.",
              srcLoc)
    {}
};

class CDT_EXPORT DuplicateVertexError : public Error
{
public:
    DuplicateVertexError(
        const VertInd v1,
        const VertInd v2,
        const SourceLocation& srcLoc)
        : Error(
              "Duplicate vertex detected: #" + CDT::to_string(v1) +
                  " is a duplicate of #" + CDT::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 CDT_EXPORT IntersectingConstraintsError : public Error
{
public:
    IntersectingConstraintsError(
        const Edge& e1,
        const Edge& e2,
        const SourceLocation& srcLoc)
        : Error(
              "Intersecting constraint edges detected: (" +
                  CDT::to_string(e1.v1()) + ", " + CDT::to_string(e1.v2()) +
                  ") intersects (" + CDT::to_string(e2.v1()) + ", " +
                  CDT::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;
};

class CDT_EXPORT InvalidEdgeSplitVertex : public Error
{
public:
    InvalidEdgeSplitVertex(
        const Edge& e1,
        const Edge& e2,
        const SourceLocation& srcLoc)
        : Error(
              "Intersection of constraint edges (" + CDT::to_string(e1.v1()) +
                  ", " + CDT::to_string(e1.v2()) + ") and (" +
                  CDT::to_string(e2.v1()) + ", " + CDT::to_string(e2.v2()) +
                  ") can not be resolved: computed split vertex is invalid",
              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;
};
 
#ifdef CDT_ENABLE_CALLBACK_HANDLER

struct CDT_EXPORT AddVertexType
{
    enum Enum
    {
        UserInput,
        FixedEdgeMidpoint,
        FixedEdgesIntersection,
    };
};

struct CDT_EXPORT TriangleChangeType
{
    enum Enum
    {
        AddedNew,         
        ModifiedExisting, 
    };
};
 
// parameter names are used for documentation purposes, even if they are un-used
// in the interface's default implementation
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"

class CDT_EXPORT ICallbackHandler
{
public:
    virtual ~ICallbackHandler()
    {}

    virtual void onAddSuperTriangle()
    {}

    virtual void onInsertVertexInsideTriangle(
        const TriInd iRepurposedTri,
        const TriInd iNewTri1,
        const TriInd iNewTri2)
    {}

    virtual void onInsertVertexOnEdge(
        const TriInd iRepurposedTri1,
        const TriInd iRepurposedTri2,
        const TriInd iNewTri1,
        const TriInd iNewTri2)
    {}

    virtual void onFlipEdge(const TriInd iT, const TriInd iTopo)
    {}

    virtual void
    onAddVertexStart(const VertInd iV, const AddVertexType::Enum vertexType)
    {}

    virtual void onAddEdgeStart(const Edge& edge)
    {}

    virtual void onReTriangulatePolygon(const std::vector<TriInd>& tris)
    {}

    virtual bool isAbortCalculation() const
    {
        return false;
    };
};
 
// parameter names are used for documentation purposes, even if they are
// un-used in the interface's default implementation
#pragma GCC diagnostic pop
 
#endif


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;
 
#ifdef CDT_ENABLE_CALLBACK_HANDLER
    void setCallbackHandler(ICallbackHandler* callbackHandler);
#endif

    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,
        const bool doHandleFixedSplitEdge = false);
    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<TriInd>& trianglesToReuse,
        std::vector<TriangulatePseudoPolygonTask>& iterations);
    void triangulatePseudoPolygonIteration(
        const std::vector<VertInd>& poly,
        unordered_map<Edge, TriInd>& outerTris,
        std::vector<TriInd>& trianglesToReuse,
        std::vector<TriangulatePseudoPolygonTask>& iterations);
    IndexSizeType findDelaunayPoint(
        const std::vector<VertInd>& poly,
        IndexSizeType iA,
        IndexSizeType iB) const;
    TriInd addTriangle(const Triangle& t);
    TriInd addTriangle();
    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);
    bool isEdgeSplitVertexValid(
        const V2d<T>& splitVert,
        TriInd iT,
        TriInd iTopo,
        VertInd iVL,
        VertInd iVR) const;
    Edge originalInputEdge(const Edge& e) const;
    const Triangle& triangleAt(TriInd iT) const;
    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, Box2d<T> box);
    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();
 
    TNearPointLocator m_nearPtLocator;
    VertInd m_nTargetVerts;
    SuperGeometryType::Enum m_superGeomType;
    VertexInsertionOrder::Enum m_vertexInsertionOrder;
    IntersectingConstraintEdges::Enum m_intersectingEdgesStrategy;
    T m_minDistToConstraintEdge;
    TriIndVec m_vertTris; 
#ifdef CDT_ENABLE_CALLBACK_HANDLER
    ICallbackHandler* m_callbackHandler; 
#endif
};

 
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)]);
    }
}

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())
        handleException(FinalizedError(CDT_SOURCE_LOCATION));
 
    const bool isFirstTime = vertices.empty();
 
    //
    // performance optimization: pre-allocate triangles and vertices
    //
    const std::size_t nNewVertices = std::distance(first, last);
    std::size_t exactCapacityTriangles = triangles.size() + 2 * nNewVertices;
    std::size_t exactCapacityVertices = vertices.size() + nNewVertices;
    if(isFirstTime) // account for adding super-triangle on the first run
    {
        exactCapacityTriangles += 1;
        exactCapacityVertices += nSuperTriangleVertices;
    }
    std::size_t capacityTriangles = exactCapacityTriangles;
    std::size_t capacityVertices = exactCapacityVertices;
    // to avoid re-allocation and unused memory
    // over-allocate by a fixed factor
    // when constraint edge intersections are resolved and vertices are many
    // because vertex is added for each intersection
    // and total number of intersections is unknown
    const VertInd overAllocationVerticesThreshold(1000);
    const T overAllocationFactor(1.1);
    const bool isOverPreAllocated =
        m_intersectingEdgesStrategy ==
            IntersectingConstraintEdges::TryResolve &&
        VertInd(nNewVertices) >= overAllocationVerticesThreshold;
    if(isOverPreAllocated)
    {
        capacityTriangles *= overAllocationFactor;
        capacityVertices *= overAllocationFactor;
    }
    triangles.reserve(capacityTriangles);
    vertices.reserve(capacityVertices);
    m_vertTris.reserve(capacityVertices);
 
    Box2d<T> box;
    if(isFirstTime)
    {
        box.envelopPoints(first, last, getX, getY);
        addSuperTriangle(box);
    }
    tryInitNearestPointLocator();
    const VertInd nExistingVerts = static_cast<VertInd>(vertices.size());
 
    for(TVertexIter it = first; it != last; ++it)
        addNewVertex(V2d<T>(getX(*it), getY(*it)), noNeighbor);
 
    switch(m_vertexInsertionOrder)
    {
    case VertexInsertionOrder::AsProvided:
        insertVertices_AsProvided(nExistingVerts);
        break;
    case VertexInsertionOrder::Auto:
        isFirstTime ? insertVertices_KDTreeBFS(nExistingVerts, box)
                    : insertVertices_Randomized(nExistingVerts);
        break;
    }
 
// make sure pre-allocation was correct
#ifdef CDT_ENABLE_CALLBACK_HANDLER
    assert(
        !m_callbackHandler || m_callbackHandler->isAbortCalculation() ||
        (vertices.size() == exactCapacityVertices));
    assert(
        !m_callbackHandler || m_callbackHandler->isAbortCalculation() ||
        (triangles.size() == exactCapacityTriangles));
#else
    assert(vertices.size() == exactCapacityVertices);
    assert(triangles.size() == exactCapacityTriangles);
#endif
}
 
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())
        handleException(FinalizedError(CDT_SOURCE_LOCATION));
 
    std::vector<TriangulatePseudoPolygonTask> tppIterations;
    EdgeVec remaining;
    for(; first != last; ++first)
    {
#ifdef CDT_ENABLE_CALLBACK_HANDLER
        if(m_callbackHandler && m_callbackHandler->isAbortCalculation())
        {
            return;
        }
#endif
        // +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);
    }
}
 
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())
        handleException(FinalizedError(CDT_SOURCE_LOCATION));
 
    tryInitNearestPointLocator();
    // state shared between different runs for performance gains
    std::vector<ConformToEdgeTask> remaining;
    for(; first != last; ++first)
    {
#ifdef CDT_ENABLE_CALLBACK_HANDLER
        if(m_callbackHandler && m_callbackHandler->isAbortCalculation())
        {
            return;
        }
#endif
        // +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);
    }
}
 
} // namespace CDT
 
#ifndef CDT_USE_AS_COMPILED_LIBRARY
#include "Triangulation.hpp"
#endif
 
#endif // header-guard