hcdt 0.1.0.2 → 0.1.0.3

raw patch · 12 files changed

+4995/−12 lines, 12 files

Files

• cpp/CDT.h +448/−0
• cpp/CDT.hpp +107/−0
• cpp/CDTUtils.h +475/−0
• cpp/CDTUtils.hpp +279/−0
• cpp/KDTree.h +399/−0
• cpp/LocatorKDTree.h +71/−0
• cpp/Triangulation.h +617/−0
• cpp/Triangulation.hpp +1558/−0
• cpp/hcdt.hpp +35/−0
• cpp/predicates.h +939/−0
• cpp/remove_at.hpp +55/−0
• hcdt.cabal +12/−12

@@ -0,0 +1,448 @@+/* 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/. */++/**+ * @file+ * Public API+ */++#ifndef CDT_lNrmUayWQaIR5fxnsg9B+#define CDT_lNrmUayWQaIR5fxnsg9B++#include "CDTUtils.h"+#include "Triangulation.h"++#include "remove_at.hpp"++#include <algorithm>+#include <cassert>+#include <cstdlib>+#include <iterator>+#include <memory>+#include <stack>+#include <vector>++/// Namespace containing triangulation functionality+namespace CDT+{++/** @defgroup API Public API+ *  Contains API for constrained and conforming Delaunay triangulations+ */+/// @{++/**+ * Type used for storing layer depths for triangles+ * @note LayerDepth should support 60K+ layers, which could be to much or+ * too little for some use cases. Feel free to re-define this typedef.+ */+typedef unsigned short LayerDepth;+typedef LayerDepth BoundaryOverlapCount;++/// Triangles by vertex index+typedef std::vector<TriIndVec> VerticesTriangles;++/** @defgroup helpers Helpers+ *  Helpers for working with CDT::Triangulation.+ */+/// @{++/**+ * Calculate triangles adjacent to vertices (triangles by vertex index)+ * @param triangles triangulation+ * @param verticesSize total number of vertices to pre-allocate the output+ * @return triangles by vertex index+ */+CDT_EXPORT VerticesTriangles+calculateTrianglesByVertex(const TriangleVec& triangles, VertInd verticesSize);++/**+ * Information about removed duplicated vertices.+ *+ * Contains mapping information and removed duplicates indices.+ * @note vertices {0,1,2,3,4} where 0 and 3 are the same will produce mapping+ *       {0,1,2,0,3} (to new vertices {0,1,2,3}) and duplicates {3}+ */+struct CDT_EXPORT DuplicatesInfo+{+    std::vector<std::size_t> mapping;    ///< vertex index mapping+    std::vector<std::size_t> duplicates; ///< duplicates' indices+};++/**+ * Find duplicates in given custom point-type range+ * @note duplicates are points with exactly same X and Y coordinates+ * @tparam TVertexIter iterator that dereferences to custom point type+ * @tparam TGetVertexCoordX function object getting x coordinate from vertex.+ * Getter signature: const TVertexIter::value_type& -> T+ * @tparam TGetVertexCoordY function object getting y coordinate from vertex.+ * Getter signature: const TVertexIter::value_type& -> T+ * @param first beginning of the range of vertices+ * @param last end of the range of vertices+ * @param getX getter of X-coordinate+ * @param getY getter of Y-coordinate+ * @returns information about vertex duplicates+ */+template <+    typename T,+    typename TVertexIter,+    typename TGetVertexCoordX,+    typename TGetVertexCoordY>+DuplicatesInfo FindDuplicates(+    TVertexIter first,+    TVertexIter last,+    TGetVertexCoordX getX,+    TGetVertexCoordY getY);++/**+ * Remove duplicates in-place from vector of custom points+ * @tparam TVertex vertex type+ * @tparam TAllocator allocator used by input vector of vertices+ * @param vertices vertices to remove duplicates from+ * @param duplicates information about duplicates+ */+template <typename TVertex, typename TAllocator>+void RemoveDuplicates(+    std::vector<TVertex, TAllocator>& vertices,+    const std::vector<std::size_t>& duplicates);++/**+ * Remove duplicated points in-place+ *+ * @tparam T type of vertex coordinates (e.g., float, double)+ * @param[in, out] vertices collection of vertices to remove duplicates from+ * @returns information about duplicated vertices that were removed.+ */+template <typename T>+CDT_EXPORT DuplicatesInfo RemoveDuplicates(std::vector<V2d<T> >& vertices);++/**+ * Remap vertex indices in edges (in-place) using given vertex-index mapping.+ * @tparam TEdgeIter iterator that dereferences to custom edge type+ * @tparam TGetEdgeVertexStart function object getting start vertex index+ * from an edge.+ * Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+ * @tparam TGetEdgeVertexEnd function object getting end vertex index from+ * an edge. Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+ * @tparam TMakeEdgeFromStartAndEnd function object that makes new edge from+ * start and end vertices+ * @param first beginning of the range of edges+ * @param last end of the range of edges+ * @param mapping vertex-index mapping+ * @param getStart getter of edge start vertex index+ * @param getEnd getter of edge end vertex index+ * @param makeEdge factory for making edge from vetices+ */+template <+    typename TEdgeIter,+    typename TGetEdgeVertexStart,+    typename TGetEdgeVertexEnd,+    typename TMakeEdgeFromStartAndEnd>+CDT_EXPORT void RemapEdges(+    TEdgeIter first,+    TEdgeIter last,+    const std::vector<std::size_t>& mapping,+    TGetEdgeVertexStart getStart,+    TGetEdgeVertexEnd getEnd,+    TMakeEdgeFromStartAndEnd makeEdge);++/**+ * Remap vertex indices in edges (in-place) using given vertex-index mapping.+ *+ * @note Mapping can be a result of RemoveDuplicates function+ * @param[in,out] edges collection of edges to remap+ * @param mapping vertex-index mapping+ */+CDT_EXPORT void+RemapEdges(std::vector<Edge>& edges, const std::vector<std::size_t>& mapping);++/**+ * Find point duplicates, remove them from vector (in-place) and remap edges+ * (in-place)+ * @note Same as a chained call of CDT::FindDuplicates, CDT::RemoveDuplicates,+ * and CDT::RemapEdges+ * @tparam T type of vertex coordinates (e.g., float, double)+ * @tparam TVertex type of vertex+ * @tparam TGetVertexCoordX function object getting x coordinate from vertex.+ * Getter signature: const TVertexIter::value_type& -> T+ * @tparam TGetVertexCoordY function object getting y coordinate from vertex.+ * Getter signature: const TVertexIter::value_type& -> T+ * @tparam TEdgeIter iterator that dereferences to custom edge type+ * @tparam TGetEdgeVertexStart function object getting start vertex index+ * from an edge.+ * Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+ * @tparam TGetEdgeVertexEnd function object getting end vertex index from+ * an edge. Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+ * @tparam TMakeEdgeFromStartAndEnd function object that makes new edge from+ * start and end vertices+ * @param[in, out] vertices vertices to remove duplicates from+ * @param[in, out] edges collection of edges connecting vertices+ * @param getX getter of X-coordinate+ * @param getY getter of Y-coordinate+ * @param edgesFirst beginning of the range of edges+ * @param edgesLast end of the range of edges+ * @param getStart getter of edge start vertex index+ * @param getEnd getter of edge end vertex index+ * @param makeEdge factory for making edge from vetices+ * @returns information about vertex duplicates+ */+template <+    typename T,+    typename TVertex,+    typename TGetVertexCoordX,+    typename TGetVertexCoordY,+    typename TVertexAllocator,+    typename TEdgeIter,+    typename TGetEdgeVertexStart,+    typename TGetEdgeVertexEnd,+    typename TMakeEdgeFromStartAndEnd>+DuplicatesInfo RemoveDuplicatesAndRemapEdges(+    std::vector<TVertex, TVertexAllocator>& vertices,+    TGetVertexCoordX getX,+    TGetVertexCoordY getY,+    TEdgeIter edgesFirst,+    TEdgeIter edgesLast,+    TGetEdgeVertexStart getStart,+    TGetEdgeVertexEnd getEnd,+    TMakeEdgeFromStartAndEnd makeEdge);++/**+ * Same as a chained call of CDT::RemoveDuplicates + CDT::RemapEdges+ *+ * @tparam T type of vertex coordinates (e.g., float, double)+ * @param[in, out] vertices collection of vertices to remove duplicates from+ * @param[in,out] edges collection of edges to remap+ */+template <typename T>+CDT_EXPORT DuplicatesInfo RemoveDuplicatesAndRemapEdges(+    std::vector<V2d<T> >& vertices,+    std::vector<Edge>& edges);++/**+ * Extract all edges of triangles+ *+ * @param triangles triangles used to extract edges+ * @return an unordered set of all edges of triangulation+ */+CDT_EXPORT EdgeUSet extractEdgesFromTriangles(const TriangleVec& triangles);++/*!+ * Converts piece->original_edges mapping to original_edge->pieces+ * @param pieceToOriginals maps pieces to original edges+ * @return mapping of original edges to pieces+ */+CDT_EXPORT unordered_map<Edge, EdgeVec>+EdgeToPiecesMapping(const unordered_map<Edge, EdgeVec>& pieceToOriginals);++/*!+ * Convert edge-to-pieces mapping into edge-to-split-vertices mapping+ * @tparam T type of vertex coordinates (e.g., float, double)+ * @param edgeToPieces edge-to-pieces mapping+ * @param vertices vertex buffer+ * @return mapping of edge-to-split-points.+ * Split points are sorted from edge's start (v1) to end (v2)+ */+template <typename T>+CDT_EXPORT unordered_map<Edge, std::vector<VertInd> > EdgeToSplitVertices(+    const unordered_map<Edge, EdgeVec>& edgeToPieces,+    const std::vector<V2d<T> >& vertices);++/// @}++/// @}++} // namespace CDT++//*****************************************************************************+// Implementations of template functionlity+//*****************************************************************************+// hash for CDT::V2d<T>+#ifdef CDT_CXX11_IS_SUPPORTED+namespace std+#else+namespace boost+#endif+{+template <typename T>+struct hash<CDT::V2d<T> >+{+    size_t operator()(const CDT::V2d<T>& xy) const+    {+#ifdef CDT_CXX11_IS_SUPPORTED+        typedef std::hash<T> Hasher;+#else+        typedef boost::hash<T> Hasher;+#endif+        return Hasher()(xy.x) ^ Hasher()(xy.y);+    }+};+} // namespace std++namespace CDT+{++//-----+// API+//-----+template <+    typename T,+    typename TVertexIter,+    typename TGetVertexCoordX,+    typename TGetVertexCoordY>+DuplicatesInfo FindDuplicates(+    TVertexIter first,+    TVertexIter last,+    TGetVertexCoordX getX,+    TGetVertexCoordY getY)+{+    typedef unordered_map<V2d<T>, std::size_t> PosToIndex;+    PosToIndex uniqueVerts;+    const std::size_t verticesSize = std::distance(first, last);+    DuplicatesInfo di = {+        std::vector<std::size_t>(verticesSize), std::vector<std::size_t>()};+    for(std::size_t iIn = 0, iOut = iIn; iIn < verticesSize; ++iIn, ++first)+    {+        typename PosToIndex::const_iterator it;+        bool isUnique;+        tie(it, isUnique) = uniqueVerts.insert(+            std::make_pair(V2d<T>::make(getX(*first), getY(*first)), iOut));+        if(isUnique)+        {+            di.mapping[iIn] = iOut++;+            continue;+        }+        di.mapping[iIn] = it->second; // found a duplicate+        di.duplicates.push_back(iIn);+    }+    return di;+}++template <typename TVertex, typename TAllocator>+void RemoveDuplicates(+    std::vector<TVertex, TAllocator>& vertices,+    const std::vector<std::size_t>& duplicates)+{+    vertices.erase(+        remove_at(+            vertices.begin(),+            vertices.end(),+            duplicates.begin(),+            duplicates.end()),+        vertices.end());+}++template <+    typename TEdgeIter,+    typename TGetEdgeVertexStart,+    typename TGetEdgeVertexEnd,+    typename TMakeEdgeFromStartAndEnd>+void RemapEdges(+    TEdgeIter first,+    const TEdgeIter last,+    const std::vector<std::size_t>& mapping,+    TGetEdgeVertexStart getStart,+    TGetEdgeVertexEnd getEnd,+    TMakeEdgeFromStartAndEnd makeEdge)+{+    for(; first != last; ++first)+    {+        *first = makeEdge(+            static_cast<VertInd>(mapping[getStart(*first)]),+            static_cast<VertInd>(mapping[getEnd(*first)]));+    }+}++template <+    typename T,+    typename TVertex,+    typename TGetVertexCoordX,+    typename TGetVertexCoordY,+    typename TVertexAllocator,+    typename TEdgeIter,+    typename TGetEdgeVertexStart,+    typename TGetEdgeVertexEnd,+    typename TMakeEdgeFromStartAndEnd>+DuplicatesInfo RemoveDuplicatesAndRemapEdges(+    std::vector<TVertex, TVertexAllocator>& vertices,+    TGetVertexCoordX getX,+    TGetVertexCoordY getY,+    const TEdgeIter edgesFirst,+    const TEdgeIter edgesLast,+    TGetEdgeVertexStart getStart,+    TGetEdgeVertexEnd getEnd,+    TMakeEdgeFromStartAndEnd makeEdge)+{+    const DuplicatesInfo di =+        FindDuplicates<T>(vertices.begin(), vertices.end(), getX, getY);+    RemoveDuplicates(vertices, di.duplicates);+    RemapEdges(edgesFirst, edgesLast, di.mapping, getStart, getEnd, makeEdge);+    return di;+}++template <typename T>+unordered_map<Edge, std::vector<VertInd> > EdgeToSplitVertices(+    const unordered_map<Edge, EdgeVec>& edgeToPieces,+    const std::vector<V2d<T> >& vertices)+{+    typedef std::pair<VertInd, T> VertCoordPair;+    struct ComparePred+    {+        bool operator()(const VertCoordPair& a, const VertCoordPair& b) const+        {+            return a.second < b.second;+        }+    } comparePred;++    unordered_map<Edge, std::vector<VertInd> > edgeToSplitVerts;+    typedef unordered_map<Edge, EdgeVec>::const_iterator It;+    for(It it = edgeToPieces.begin(); it != edgeToPieces.end(); ++it)+    {+        const Edge& e = it->first;+        const T dX = vertices[e.v2()].x - vertices[e.v1()].x;+        const T dY = vertices[e.v2()].y - vertices[e.v1()].y;+        const bool isX = std::abs(dX) >= std::abs(dY); // X-coord longer+        const bool isAscending =+            isX ? dX >= 0 : dY >= 0; // Longer coordinate ascends+        const EdgeVec& pieces = it->second;+        std::vector<VertCoordPair> splitVerts;+        // size is:  2[ends] + (pieces - 1)[split vertices] = pieces + 1+        splitVerts.reserve(pieces.size() + 1);+        typedef EdgeVec::const_iterator EIt;+        for(EIt it = pieces.begin(); it != pieces.end(); ++it)+        {+            const array<VertInd, 2> vv = {it->v1(), it->v2()};+            typedef array<VertInd, 2>::const_iterator VIt;+            for(VIt v = vv.begin(); v != vv.end(); ++v)+            {+                const T c = isX ? vertices[*v].x : vertices[*v].y;+                splitVerts.push_back(std::make_pair(*v, isAscending ? c : -c));+            }+        }+        // sort by longest coordinate+        std::sort(splitVerts.begin(), splitVerts.end(), comparePred);+        // remove duplicates+        splitVerts.erase(+            std::unique(splitVerts.begin(), splitVerts.end()),+            splitVerts.end());+        assert(splitVerts.size() > 2); // 2 end points with split vertices+        std::pair<Edge, std::vector<VertInd> > val =+            std::make_pair(e, std::vector<VertInd>());+        val.second.reserve(splitVerts.size());+        typedef typename std::vector<VertCoordPair>::const_iterator SEIt;+        for(SEIt it = splitVerts.begin() + 1; it != splitVerts.end() - 1; ++it)+        {+            val.second.push_back(it->first);+        }+        edgeToSplitVerts.insert(val);+    }+    return edgeToSplitVerts;+}++} // namespace CDT++#ifndef CDT_USE_AS_COMPILED_LIBRARY+#include "CDT.hpp"+#endif++#endif // header-guard

@@ -0,0 +1,107 @@+/* 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/. */++/**+ * @file+ * Public API - implementation+ */++#include "CDT.h"++#include <algorithm>+#include <deque>+#include <limits>+#include <stdexcept>++namespace CDT+{++CDT_INLINE_IF_HEADER_ONLY VerticesTriangles calculateTrianglesByVertex(+    const TriangleVec& triangles,+    const VertInd verticesSize)+{+    VerticesTriangles vertTris(verticesSize);+    for(TriInd iT = 0; iT < triangles.size(); ++iT)+    {+        const VerticesArr3& vv = triangles[iT].vertices;+        for(VerticesArr3::const_iterator v = vv.begin(); v != vv.end(); ++v)+        {+            vertTris[*v].push_back(iT);+        }+    }+    return vertTris;+}++template <typename T>+DuplicatesInfo RemoveDuplicates(std::vector<V2d<T> >& vertices)+{+    const DuplicatesInfo di = FindDuplicates<T>(+        vertices.begin(), vertices.end(), getX_V2d<T>, getY_V2d<T>);+    RemoveDuplicates(vertices, di.duplicates);+    return di;+}++CDT_INLINE_IF_HEADER_ONLY void+RemapEdges(std::vector<Edge>& edges, const std::vector<std::size_t>& mapping)+{+    RemapEdges(+        edges.begin(),+        edges.end(),+        mapping,+        edge_get_v1,+        edge_get_v2,+        edge_make);+}++template <typename T>+DuplicatesInfo RemoveDuplicatesAndRemapEdges(+    std::vector<V2d<T> >& vertices,+    std::vector<Edge>& edges)+{+    return RemoveDuplicatesAndRemapEdges<T>(+        vertices,+        getX_V2d<T>,+        getY_V2d<T>,+        edges.begin(),+        edges.end(),+        edge_get_v1,+        edge_get_v2,+        edge_make);+}++CDT_INLINE_IF_HEADER_ONLY EdgeUSet+extractEdgesFromTriangles(const TriangleVec& triangles)+{+    EdgeUSet edges;+    typedef TriangleVec::const_iterator CIt;+    for(CIt t = triangles.begin(); t != triangles.end(); ++t)+    {+        edges.insert(Edge(VertInd(t->vertices[0]), VertInd(t->vertices[1])));+        edges.insert(Edge(VertInd(t->vertices[1]), VertInd(t->vertices[2])));+        edges.insert(Edge(VertInd(t->vertices[2]), VertInd(t->vertices[0])));+    }+    return edges;+}++CDT_INLINE_IF_HEADER_ONLY unordered_map<Edge, EdgeVec>+EdgeToPiecesMapping(const unordered_map<Edge, EdgeVec>& pieceToOriginals)+{+    unordered_map<Edge, EdgeVec> originalToPieces;+    typedef unordered_map<Edge, EdgeVec>::const_iterator Cit;+    for(Cit ptoIt = pieceToOriginals.begin(); ptoIt != pieceToOriginals.end();+        ++ptoIt)+    {+        const Edge piece = ptoIt->first;+        const EdgeVec& originals = ptoIt->second;+        for(EdgeVec::const_iterator origIt = originals.begin();+            origIt != originals.end();+            ++origIt)+        {+            originalToPieces[*origIt].push_back(piece);+        }+    }+    return originalToPieces;+}++} // namespace CDT

@@ -0,0 +1,475 @@+/* 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/. */++/**+ * @file+ * Utilities and helpers+ */++#ifndef CDT_obwOaxOTdAWcLNTlNnaq+#define CDT_obwOaxOTdAWcLNTlNnaq++#ifdef CDT_DONT_USE_BOOST_RTREE+// CDT_DONT_USE_BOOST_RTREE was replaced with CDT_USE_BOOST+typedef char CDT_DONT_USE_BOOST_RTREE__was__replaced__with__CDT_USE_BOOST[-1];+#endif++// #define CDT_USE_STRONG_TYPING // strong type checks on indices++// check if c++11 is supported+#if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)+#define CDT_CXX11_IS_SUPPORTED+#elif !defined(__cplusplus) && !defined(_MSC_VER)+typedef char couldnt_parse_cxx_standard[-1]; ///< Error: couldn't parse standard+#endif++// Functions defined outside the class need to be 'inline'+// if CDT is configured to be used as header-only library:+// single-definition rule is violated otherwise+#ifdef CDT_USE_AS_COMPILED_LIBRARY+#define CDT_INLINE_IF_HEADER_ONLY+#include "cdt_export.h" // automatically generated by CMake+#else+/**+ * Macro for inlining non-template functions when in header-only mode to+ * avoid multiple declaration errors.+ */+#define CDT_INLINE_IF_HEADER_ONLY inline+/// Export not needed in header-only mode+#define CDT_EXPORT+#endif++#include <cassert>+#include <cmath>+#include <limits>+#include <vector>++#ifdef CDT_USE_BOOST+#include <boost/container/flat_set.hpp>+#endif++// use fall-backs for c++11 features+#ifdef CDT_CXX11_IS_SUPPORTED++#include <array>+#include <functional>+#include <random>+#include <tuple>+#include <unordered_map>+#include <unordered_set>+namespace CDT+{+using std::array;+using std::get;+using std::make_tuple;+using std::mt19937;+using std::tie;+using std::tuple;+using std::unordered_map;+using std::unordered_set;+} // namespace CDT++#else+#include <boost/array.hpp>+#include <boost/functional/hash.hpp>+#include <boost/random.hpp>+#include <boost/tuple/tuple.hpp>+#include <boost/unordered_map.hpp>+#include <boost/unordered_set.hpp>+namespace CDT+{+using boost::array;+using boost::get;+using boost::make_tuple;+using boost::tie;+using boost::tuple;+using boost::unordered_map;+using boost::unordered_set;+using boost::random::mt19937;+} // namespace CDT+#endif++namespace CDT+{++/// 2D vector+template <typename T>+struct CDT_EXPORT V2d+{+    T x; ///< X-coordinate+    T y; ///< Y-coordinate++    /// Create vector from X and Y coordinates+    static V2d make(T x, T y);+};++/// X- coordinate getter for V2d+template <typename T>+const T& getX_V2d(const V2d<T>& v)+{+    return v.x;+}++/// Y-coordinate getter for V2d+template <typename T>+const T& getY_V2d(const V2d<T>& v)+{+    return v.y;+}++/// If two 2D vectors are exactly equal+template <typename T>+bool operator==(const CDT::V2d<T>& lhs, const CDT::V2d<T>& rhs)+{+    return lhs.x == rhs.x && lhs.y == rhs.y;+}++#ifdef CDT_USE_64_BIT_INDEX_TYPE+typedef unsigned long long IndexSizeType;+#else+typedef unsigned int IndexSizeType;+#endif++#ifndef CDT_USE_STRONG_TYPING+/// Index in triangle+typedef unsigned char Index;+/// Vertex index+typedef IndexSizeType VertInd;+/// Triangle index+typedef IndexSizeType TriInd;+#else+/// Index in triangle+BOOST_STRONG_TYPEDEF(unsigned char, Index);+/// Vertex index+BOOST_STRONG_TYPEDEF(IndexSizeType, VertInd);+/// Triangle index+BOOST_STRONG_TYPEDEF(IndexSizeType, TriInd);+#endif++/// Constant representing no valid neighbor for a triangle+const static TriInd noNeighbor(std::numeric_limits<TriInd>::max());+/// Constant representing no valid vertex for a triangle+const static VertInd noVertex(std::numeric_limits<VertInd>::max());++typedef std::vector<TriInd> TriIndVec;  ///< Vector of triangle indices+typedef array<VertInd, 3> VerticesArr3; ///< array of three vertex indices+typedef array<TriInd, 3> NeighborsArr3; ///< array of three neighbors++/// 2D bounding box+template <typename T>+struct CDT_EXPORT Box2d+{+    V2d<T> min; ///< min box corner+    V2d<T> max; ///< max box corner++    /// Envelop box around a point+    void envelopPoint(const V2d<T>& p)+    {+        envelopPoint(p.x, p.y);+    }+    /// Envelop box around a point with given coordinates+    void envelopPoint(const T x, const T y)+    {+        min.x = std::min(x, min.x);+        max.x = std::max(x, max.x);+        min.y = std::min(y, min.y);+        max.y = std::max(y, max.y);+    }+};++/// Bounding box of a collection of custom 2D points given coordinate getters+template <+    typename T,+    typename TVertexIter,+    typename TGetVertexCoordX,+    typename TGetVertexCoordY>+Box2d<T> envelopBox(+    TVertexIter first,+    TVertexIter last,+    TGetVertexCoordX getX,+    TGetVertexCoordY getY)+{+    const T max = std::numeric_limits<T>::max();+    Box2d<T> box = {{max, max}, {-max, -max}};+    for(; first != last; ++first)+    {+        box.envelopPoint(getX(*first), getY(*first));+    }+    return box;+}++/// Bounding box of a collection of 2D points+template <typename T>+CDT_EXPORT Box2d<T> envelopBox(const std::vector<V2d<T> >& vertices);++/// Edge connecting two vertices: vertex with smaller index is always first+/// \note: hash Edge is specialized at the bottom+struct CDT_EXPORT Edge+{+    /// Constructor+    Edge(VertInd iV1, VertInd iV2);+    /// Equals operator+    bool operator==(const Edge& other) const;+    /// Not-equals operator+    bool operator!=(const Edge& other) const;+    /// V1 getter+    VertInd v1() const;+    /// V2 getter+    VertInd v2() const;+    /// Edges' vertices+    const std::pair<VertInd, VertInd>& verts() const;++private:+    std::pair<VertInd, VertInd> m_vertices;+};++/// Get edge first vertex+inline VertInd edge_get_v1(const Edge& e)+{+    return e.v1();+}++/// Get edge second vertex+inline VertInd edge_get_v2(const Edge& e)+{+    return e.v2();+}++/// Get edge second vertex+inline Edge edge_make(VertInd iV1, VertInd iV2)+{+    return Edge(iV1, iV2);+}++typedef std::vector<Edge> EdgeVec;                ///< Vector of edges+typedef unordered_set<Edge> EdgeUSet;             ///< Hash table of edges+typedef unordered_set<TriInd> TriIndUSet;         ///< Hash table of triangles+typedef unordered_map<TriInd, TriInd> TriIndUMap; ///< Triangle hash map+#ifdef CDT_USE_BOOST+/// Flat hash table of triangles+typedef boost::container::flat_set<TriInd> TriIndFlatUSet;+#endif++/// Triangulation triangle (CCW winding)+/* Counter-clockwise winding:+       v3+       /\+    n3/  \n2+     /____\+   v1  n1  v2                 */+struct CDT_EXPORT Triangle+{+    VerticesArr3 vertices;   ///< triangle's three vertices+    NeighborsArr3 neighbors; ///< triangle's three neighbors++    /**+     * Factory method+     * @note needed for c++03 compatibility (no uniform initialization+     * available)+     */+    static Triangle+    make(const array<VertInd, 3>& vertices, const array<TriInd, 3>& neighbors)+    {+        Triangle t = {vertices, neighbors};+        return t;+    }+};++typedef std::vector<Triangle> TriangleVec; ///< Vector of triangles++/// Advance vertex or neighbor index counter-clockwise+CDT_EXPORT Index ccw(Index i);++/// Advance vertex or neighbor index clockwise+CDT_EXPORT Index cw(Index i);++/// Location of point on a triangle+struct CDT_EXPORT PtTriLocation+{+    /// Enum+    enum Enum+    {+        Inside,+        Outside,+        OnEdge1,+        OnEdge2,+        OnEdge3,+    };+};++/// Check if location is classified as on any of three edges+CDT_EXPORT bool isOnEdge(PtTriLocation::Enum location);++/// Neighbor index from a on-edge location+/// \note Call only if located on the edge!+CDT_EXPORT Index edgeNeighbor(PtTriLocation::Enum location);++/// Relative location of point to a line+struct CDT_EXPORT PtLineLocation+{+    /// Enum+    enum Enum+    {+        Left,+        Right,+        OnLine,+    };+};++/// Orient p against line v1-v2 2D: robust geometric predicate+template <typename T>+CDT_EXPORT T orient2D(const V2d<T>& p, const V2d<T>& v1, const V2d<T>& v2);++/// Check if point lies to the left of, to the right of, or on a line+template <typename T>+CDT_EXPORT PtLineLocation::Enum locatePointLine(+    const V2d<T>& p,+    const V2d<T>& v1,+    const V2d<T>& v2,+    T orientationTolerance = T(0));++/// Classify value of orient2d predicate+template <typename T>+CDT_EXPORT PtLineLocation::Enum+classifyOrientation(T orientation, T orientationTolerance = T(0));++/// Check if point a lies inside of, outside of, or on an edge of a triangle+template <typename T>+CDT_EXPORT PtTriLocation::Enum locatePointTriangle(+    const V2d<T>& p,+    const V2d<T>& v1,+    const V2d<T>& v2,+    const V2d<T>& v3);++/// Opposed neighbor index from vertex index+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index opoNbr(Index vertIndex);++/// Opposed vertex index from neighbor index+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index opoVrt(Index neighborIndex);++/// Index of triangle's neighbor opposed to a vertex+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index+opposedTriangleInd(const Triangle& tri, VertInd iVert);++/// Index of triangle's neighbor opposed to an edge+CDT_INLINE_IF_HEADER_ONLY Index+opposedTriangleInd(const Triangle& tri, VertInd iVedge1, VertInd iVedge2);++/// Index of triangle's vertex opposed to a triangle+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index+opposedVertexInd(const Triangle& tri, TriInd iTopo);++/// If triangle has a given neighbor return neighbor-index, throw otherwise+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index+neighborInd(const Triangle& tri, TriInd iTnbr);++/// If triangle has a given vertex return vertex-index, throw otherwise+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY Index+vertexInd(const Triangle& tri, VertInd iV);++/// Given triangle and a vertex find opposed triangle+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY TriInd+opposedTriangle(const Triangle& tri, VertInd iVert);++/// Given two triangles, return vertex of first triangle opposed to the second+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY VertInd+opposedVertex(const Triangle& tri, TriInd iTopo);++/// Test if point lies in a circumscribed circle of a triangle+template <typename T>+CDT_EXPORT bool isInCircumcircle(+    const V2d<T>& p,+    const V2d<T>& v1,+    const V2d<T>& v2,+    const V2d<T>& v3);++/// Test if two vertices share at least one common triangle+CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY bool+verticesShareEdge(const TriIndVec& aTris, const TriIndVec& bTris);++/// Distance between two 2D points+template <typename T>+CDT_EXPORT T distance(const V2d<T>& a, const V2d<T>& b);++/// Squared distance between two 2D points+template <typename T>+CDT_EXPORT T distanceSquared(const V2d<T>& a, const V2d<T>& b);++} // namespace CDT++#ifndef CDT_USE_AS_COMPILED_LIBRARY+#include "CDTUtils.hpp"+#endif++//*****************************************************************************+// Specialize hash functions+//*****************************************************************************+#ifdef CDT_CXX11_IS_SUPPORTED+namespace std+#else+namespace boost+#endif+{++#ifdef CDT_USE_STRONG_TYPING++/// Vertex index hasher+template <>+struct hash<CDT::VertInd>+{+    /// Hash operator+    std::size_t operator()(const CDT::VertInd& vi) const+    {+        return std::hash<std::size_t>()(vi.t);+    }+};++/// Triangle index hasher+template <>+struct hash<CDT::TriInd>+{+    /// Hash operator+    std::size_t operator()(const CDT::TriInd& vi) const+    {+        return std::hash<std::size_t>()(vi.t);+    }+};++#endif // CDT_USE_STRONG_TYPING++/// Edge hasher+template <>+struct hash<CDT::Edge>+{+    /// Hash operator+    std::size_t operator()(const CDT::Edge& e) const+    {+        return hashEdge(e);+    }++private:+    static void hashCombine(std::size_t& seed, const CDT::VertInd& key)+    {+#ifdef CDT_CXX11_IS_SUPPORTED+        typedef std::hash<CDT::VertInd> Hasher;+#else+        typedef boost::hash<CDT::VertInd> Hasher;+#endif+        seed ^= Hasher()(key) + 0x9e3779b9 + (seed << 6) + (seed >> 2);+    }+    static std::size_t hashEdge(const CDT::Edge& e)+    {+        const std::pair<CDT::VertInd, CDT::VertInd>& vv = e.verts();+        std::size_t seed1(0);+        hashCombine(seed1, vv.first);+        hashCombine(seed1, vv.second);+        std::size_t seed2(0);+        hashCombine(seed2, vv.second);+        hashCombine(seed2, vv.first);+        return std::min(seed1, seed2);+    }+};+} // namespace std/boost++#endif // header guard

@@ -0,0 +1,279 @@+/* 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/. */++/**+ * @file+ * Utilities and helpers - implementation+ */++#include "CDTUtils.h"++#include "predicates.h" // robust predicates: orient, in-circle++#include <stdexcept>++namespace CDT+{++//*****************************************************************************+// V2d+//*****************************************************************************+template <typename T>+V2d<T> V2d<T>::make(const T x, const T y)+{+    V2d<T> out = {x, y};+    return out;+}++//*****************************************************************************+// Box2d+//*****************************************************************************+template <typename T>+Box2d<T> envelopBox(const std::vector<V2d<T> >& vertices)+{+    return envelopBox<T>(+        vertices.begin(), vertices.end(), getX_V2d<T>, getY_V2d<T>);+}++//*****************************************************************************+// Edge+//*****************************************************************************+CDT_INLINE_IF_HEADER_ONLY Edge::Edge(VertInd iV1, VertInd iV2)+    : m_vertices(+          iV1 < iV2 ? std::make_pair(iV1, iV2) : std::make_pair(iV2, iV1))+{}++CDT_INLINE_IF_HEADER_ONLY bool Edge::operator==(const Edge& other) const+{+    return m_vertices == other.m_vertices;+}++CDT_INLINE_IF_HEADER_ONLY bool Edge::operator!=(const Edge& other) const+{+    return !(this->operator==(other));+}++CDT_INLINE_IF_HEADER_ONLY VertInd Edge::v1() const+{+    return m_vertices.first;+}++CDT_INLINE_IF_HEADER_ONLY VertInd Edge::v2() const+{+    return m_vertices.second;+}++CDT_INLINE_IF_HEADER_ONLY const std::pair<VertInd, VertInd>& Edge::verts() const+{+    return m_vertices;+}++//*****************************************************************************+// Utility functions+//*****************************************************************************+CDT_INLINE_IF_HEADER_ONLY Index ccw(Index i)+{+    return Index((i + 1) % 3);+}++CDT_INLINE_IF_HEADER_ONLY Index cw(Index i)+{+    return Index((i + 2) % 3);+}++CDT_INLINE_IF_HEADER_ONLY bool isOnEdge(const PtTriLocation::Enum location)+{+    return location > PtTriLocation::Outside;+}++CDT_INLINE_IF_HEADER_ONLY Index edgeNeighbor(const PtTriLocation::Enum location)+{+    assert(location >= PtTriLocation::OnEdge1);+    return static_cast<Index>(location - PtTriLocation::OnEdge1);+}++template <typename T>+T orient2D(const V2d<T>& p, const V2d<T>& v1, const V2d<T>& v2)+{+    return predicates::adaptive::orient2d(v1.x, v1.y, v2.x, v2.y, p.x, p.y);+}++template <typename T>+PtLineLocation::Enum locatePointLine(+    const V2d<T>& p,+    const V2d<T>& v1,+    const V2d<T>& v2,+    const T orientationTolerance)+{+    return classifyOrientation(orient2D(p, v1, v2), orientationTolerance);+}++template <typename T>+PtLineLocation::Enum+classifyOrientation(const T orientation, const T orientationTolerance)+{+    if(orientation < -orientationTolerance)+        return PtLineLocation::Right;+    if(orientation > orientationTolerance)+        return PtLineLocation::Left;+    return PtLineLocation::OnLine;+}++template <typename T>+PtTriLocation::Enum locatePointTriangle(+    const V2d<T>& p,+    const V2d<T>& v1,+    const V2d<T>& v2,+    const V2d<T>& v3)+{+    using namespace predicates::adaptive;+    PtTriLocation::Enum result = PtTriLocation::Inside;+    PtLineLocation::Enum edgeCheck = locatePointLine(p, v1, v2);+    if(edgeCheck == PtLineLocation::Right)+        return PtTriLocation::Outside;+    if(edgeCheck == PtLineLocation::OnLine)+        result = PtTriLocation::OnEdge1;+    edgeCheck = locatePointLine(p, v2, v3);+    if(edgeCheck == PtLineLocation::Right)+        return PtTriLocation::Outside;+    if(edgeCheck == PtLineLocation::OnLine)+        result = PtTriLocation::OnEdge2;+    edgeCheck = locatePointLine(p, v3, v1);+    if(edgeCheck == PtLineLocation::Right)+        return PtTriLocation::Outside;+    if(edgeCheck == PtLineLocation::OnLine)+        result = PtTriLocation::OnEdge3;+    return result;+}++CDT_INLINE_IF_HEADER_ONLY Index opoNbr(const Index vertIndex)+{+    if(vertIndex == Index(0))+        return Index(1);+    if(vertIndex == Index(1))+        return Index(2);+    if(vertIndex == Index(2))+        return Index(0);+    throw std::runtime_error("Invalid vertex index");+}++CDT_INLINE_IF_HEADER_ONLY Index opoVrt(const Index neighborIndex)+{+    if(neighborIndex == Index(0))+        return Index(2);+    if(neighborIndex == Index(1))+        return Index(0);+    if(neighborIndex == Index(2))+        return Index(1);+    throw std::runtime_error("Invalid neighbor index");+}++CDT_INLINE_IF_HEADER_ONLY Index+opposedTriangleInd(const Triangle& tri, const VertInd iVert)+{+    for(Index vi = Index(0); vi < Index(3); ++vi)+        if(iVert == tri.vertices[vi])+            return opoNbr(vi);+    throw std::runtime_error("Could not find opposed triangle index");+}++CDT_INLINE_IF_HEADER_ONLY Index opposedTriangleInd(+    const Triangle& tri,+    const VertInd iVedge1,+    const VertInd iVedge2)+{+    for(Index vi = Index(0); vi < Index(3); ++vi)+    {+        const VertInd iVert = tri.vertices[vi];+        if(iVert != iVedge1 && iVert != iVedge2)+            return opoNbr(vi);+    }+    throw std::runtime_error("Could not find opposed-to-edge triangle index");+}++CDT_INLINE_IF_HEADER_ONLY Index+opposedVertexInd(const Triangle& tri, const TriInd iTopo)+{+    for(Index ni = Index(0); ni < Index(3); ++ni)+        if(iTopo == tri.neighbors[ni])+            return opoVrt(ni);+    throw std::runtime_error("Could not find opposed vertex index");+}++CDT_INLINE_IF_HEADER_ONLY Index+neighborInd(const Triangle& tri, const TriInd iTnbr)+{+    for(Index ni = Index(0); ni < Index(3); ++ni)+        if(iTnbr == tri.neighbors[ni])+            return ni;+    throw std::runtime_error("Could not find neighbor triangle index");+}++CDT_INLINE_IF_HEADER_ONLY Index vertexInd(const Triangle& tri, const VertInd iV)+{+    for(Index i = Index(0); i < Index(3); ++i)+        if(iV == tri.vertices[i])+            return i;+    throw std::runtime_error("Could not find vertex index in triangle");+}++CDT_INLINE_IF_HEADER_ONLY TriInd+opposedTriangle(const Triangle& tri, const VertInd iVert)+{+    return tri.neighbors[opposedTriangleInd(tri, iVert)];+}++CDT_INLINE_IF_HEADER_ONLY VertInd+opposedVertex(const Triangle& tri, const TriInd iTopo)+{+    return tri.vertices[opposedVertexInd(tri, iTopo)];+}++template <typename T>+bool isInCircumcircle(+    const V2d<T>& p,+    const V2d<T>& v1,+    const V2d<T>& v2,+    const V2d<T>& v3)+{+    using namespace predicates::adaptive;+    return incircle(v1.x, v1.y, v2.x, v2.y, v3.x, v3.y, p.x, p.y) > T(0);+}++CDT_INLINE_IF_HEADER_ONLY+bool verticesShareEdge(const TriIndVec& aTris, const TriIndVec& bTris)+{+    for(TriIndVec::const_iterator it = aTris.begin(); it != aTris.end(); ++it)+        if(std::find(bTris.begin(), bTris.end(), *it) != bTris.end())+            return true;+    return false;+}++template <typename T>+T distanceSquared(const T ax, const T ay, const T bx, const T by)+{+    const T dx = bx - ax;+    const T dy = by - ay;+    return dx * dx + dy * dy;+}++template <typename T>+T distance(const T ax, const T ay, const T bx, const T by)+{+    return std::sqrt(distanceSquared(ax, ay, bx, by));+}++template <typename T>+T distance(const V2d<T>& a, const V2d<T>& b)+{+    return distance(a.x, a.y, b.x, b.y);+}++template <typename T>+T distanceSquared(const V2d<T>& a, const V2d<T>& b)+{+    return distanceSquared(a.x, a.y, b.x, b.y);+}++} // namespace CDT

@@ -0,0 +1,399 @@+/// 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/.+/// Contribution of original implementation:+/// Andre Fecteau <andre.fecteau1@gmail.com>++#ifndef KDTREE_KDTREE_H+#define KDTREE_KDTREE_H++#include "CDTUtils.h"++#include <cassert>+#include <limits>++namespace KDTree+{++struct NodeSplitDirection+{+    enum Enum+    {+        X,+        Y,+    };+};++/// Simple tree structure with alternating half splitting nodes+/// @details Simple tree structure+///          - Tree to incrementally add points to the structure.+///          - Get the nearest point to a given input.+///          - Does not check for duplicates, expect unique points.+/// @tparam TCoordType type used for storing point coordinate.+/// @tparam NumVerticesInLeaf The number of points per leaf.+/// @tparam InitialStackDepth initial size of stack depth for nearest query.+/// Should be at least 1.+/// @tparam StackDepthIncrement increment of stack depth for nearest query when+/// stack depth is reached.+template <+    typename TCoordType,+    size_t NumVerticesInLeaf,+    size_t InitialStackDepth,+    size_t StackDepthIncrement>+class KDTree+{+public:+    typedef TCoordType coord_type;+    typedef CDT::V2d<coord_type> point_type;+    typedef CDT::VertInd point_index;+    typedef std::pair<point_type, point_index> value_type;+    typedef std::vector<point_index> point_data_vec;+    typedef point_data_vec::const_iterator pd_cit;+    typedef CDT::VertInd node_index;+    typedef CDT::array<node_index, 2> children_type;++    /// Stores kd-tree node data+    struct Node+    {+        children_type children; ///< two children if not leaf; {0,0} if leaf+        point_data_vec data;    ///< points' data if leaf+        /// Create empty leaf+        Node()+        {+            setChildren(0, 0);+            data.reserve(NumVerticesInLeaf);+        }+        /// Children setter for convenience+        void setChildren(const node_index c1, const node_index c2)+        {+            children[0] = c1;+            children[1] = c2;+        }+        /// Check if node is a leaf (has no valid children)+        bool isLeaf() const+        {+            return children[0] == children[1];+        }+    };++    /// Default constructor+    KDTree()+        : m_rootDir(NodeSplitDirection::X)+        , m_min(point_type::make(+              -std::numeric_limits<coord_type>::max(),+              -std::numeric_limits<coord_type>::max()))+        , m_max(point_type::make(+              std::numeric_limits<coord_type>::max(),+              std::numeric_limits<coord_type>::max()))+        , m_isRootBoxInitialized(false)+        , m_tasksStack(InitialStackDepth, NearestTask())+    {+        m_root = addNewNode();+    }++    /// Constructor with bounding box known in advance+    KDTree(const point_type& min, const point_type& max)+        : m_rootDir(NodeSplitDirection::X)+        , m_min(min)+        , m_max(max)+        , m_isRootBoxInitialized(true)+        , m_tasksStack(InitialStackDepth, NearestTask())+    {+        m_root = addNewNode();+    }++    /// Insert a point into kd-tree+    /// @note external point-buffer is used to reduce kd-tree's memory footprint+    /// @param iPoint index of point in external point-buffer+    /// @param points external point-buffer+    void+    insert(const point_index& iPoint, const std::vector<point_type>& points)+    {+        // if point is outside root, extend tree by adding new roots+        const point_type& pos = points[iPoint];+        while(!isInsideBox(pos, m_min, m_max))+        {+            extendTree(pos);+        }+        // now insert the point into the tree+        node_index node = m_root;+        point_type min = m_min;+        point_type max = m_max;+        NodeSplitDirection::Enum dir = m_rootDir;++        // below: initialized only to suppress warnings+        NodeSplitDirection::Enum newDir(NodeSplitDirection::X);+        coord_type mid(0);+        point_type newMin, newMax;+        while(true)+        {+            if(m_nodes[node].isLeaf())+            {+                // add point if capacity is not reached+                point_data_vec& pd = m_nodes[node].data;+                if(pd.size() < NumVerticesInLeaf)+                {+                    pd.push_back(iPoint);+                    return;+                }+                // initialize bbox first time the root capacity is reached+                if(!m_isRootBoxInitialized)+                {+                    initializeRootBox(points);+                    min = m_min;+                    max = m_max;+                }+                // split a full leaf node+                calcSplitInfo(min, max, dir, mid, newDir, newMin, newMax);+                const node_index c1 = addNewNode(), c2 = addNewNode();+                Node& n = m_nodes[node];+                n.setChildren(c1, c2);+                point_data_vec& c1data = m_nodes[c1].data;+                point_data_vec& c2data = m_nodes[c2].data;+                // move node's points to children+                for(pd_cit it = n.data.begin(); it != n.data.end(); ++it)+                {+                    whichChild(points[*it], mid, dir) == 0+                        ? c1data.push_back(*it)+                        : c2data.push_back(*it);+                }+                n.data = point_data_vec();+            }+            else+            {+                calcSplitInfo(min, max, dir, mid, newDir, newMin, newMax);+            }+            // add the point to a child+            const std::size_t iChild = whichChild(points[iPoint], mid, dir);+            iChild == 0 ? max = newMax : min = newMin;+            node = m_nodes[node].children[iChild];+            dir = newDir;+        }+    }++    /// Query kd-tree for a nearest neighbor point+    /// @note external point-buffer is used to reduce kd-tree's memory footprint+    /// @param point query point position+    /// @param points external point-buffer+    value_type nearest(+        const point_type& point,+        const std::vector<point_type>& points) const+    {+        value_type out;+        int iTask = -1;+        coord_type minDistSq = std::numeric_limits<coord_type>::max();+        m_tasksStack[++iTask] =+            NearestTask(m_root, m_min, m_max, m_rootDir, minDistSq);+        while(iTask != -1)+        {+            const NearestTask t = m_tasksStack[iTask--];+            if(t.distSq > minDistSq)+                continue;+            const Node& n = m_nodes[t.node];+            if(n.isLeaf())+            {+                for(pd_cit it = n.data.begin(); it != n.data.end(); ++it)+                {+                    const point_type& p = points[*it];+                    const coord_type distSq = CDT::distanceSquared(point, p);+                    if(distSq < minDistSq)+                    {+                        minDistSq = distSq;+                        out.first = p;+                        out.second = *it;+                    }+                }+            }+            else+            {+                coord_type mid(0);+                NodeSplitDirection::Enum newDir;+                point_type newMin, newMax;+                calcSplitInfo(t.min, t.max, t.dir, mid, newDir, newMin, newMax);++                const coord_type distToMid = t.dir == NodeSplitDirection::X+                                                 ? (point.x - mid)+                                                 : (point.y - mid);+                const coord_type toMidSq = distToMid * distToMid;++                const std::size_t iChild = whichChild(point, mid, t.dir);+                if(iTask + 2 >= static_cast<int>(m_tasksStack.size()))+                {+                    m_tasksStack.resize(+                        m_tasksStack.size() + StackDepthIncrement);+                }+                // node containing point should end up on top of the stack+                if(iChild == 0)+                {+                    m_tasksStack[++iTask] = NearestTask(+                        n.children[1], newMin, t.max, newDir, toMidSq);+                    m_tasksStack[++iTask] = NearestTask(+                        n.children[0], t.min, newMax, newDir, toMidSq);+                }+                else+                {+                    m_tasksStack[++iTask] = NearestTask(+                        n.children[0], t.min, newMax, newDir, toMidSq);+                    m_tasksStack[++iTask] = NearestTask(+                        n.children[1], newMin, t.max, newDir, toMidSq);+                }+            }+        }+        return out;+    }++private:+    /// Add a new node and return it's index in nodes buffer+    node_index addNewNode()+    {+        const node_index newNodeIndex = static_cast<node_index>(m_nodes.size());+        m_nodes.push_back(Node());+        return newNodeIndex;+    }++    /// Test which child point belongs to after the split+    /// @returns 0 if first child, 1 if second child+    std::size_t whichChild(+        const point_type& point,+        const coord_type& split,+        const NodeSplitDirection::Enum dir) const+    {+        return static_cast<size_t>(+            dir == NodeSplitDirection::X ? point.x > split : point.y > split);+    }++    /// Calculate split location, direction, and children boxes+    static void calcSplitInfo(+        const point_type& min,+        const point_type& max,+        const NodeSplitDirection::Enum dir,+        coord_type& midOut,+        NodeSplitDirection::Enum& newDirOut,+        point_type& newMinOut,+        point_type& newMaxOut)+    {+        newMaxOut = max;+        newMinOut = min;+        switch(dir)+        {+        case NodeSplitDirection::X:+            midOut = (min.x + max.x) / coord_type(2);+            newDirOut = NodeSplitDirection::Y;+            newMinOut.x = midOut;+            newMaxOut.x = midOut;+            return;+        case NodeSplitDirection::Y:+            midOut = (min.y + max.y) / coord_type(2);+            newDirOut = NodeSplitDirection::X;+            newMinOut.y = midOut;+            newMaxOut.y = midOut;+            return;+        }+    }++    /// Test if point is inside a box+    static bool isInsideBox(+        const point_type& p,+        const point_type& min,+        const point_type& max)+    {+        return p.x >= min.x && p.x <= max.x && p.y >= min.y && p.y <= max.y;+    }++    /// Extend a tree by creating new root with old root and a new node as+    /// children+    void extendTree(const point_type& point)+    {+        const node_index newRoot = addNewNode();+        const node_index newLeaf = addNewNode();+        switch(m_rootDir)+        {+        case NodeSplitDirection::X:+            m_rootDir = NodeSplitDirection::Y;+            point.y < m_min.y ? m_nodes[newRoot].setChildren(newLeaf, m_root)+                              : m_nodes[newRoot].setChildren(m_root, newLeaf);+            if(point.y < m_min.y)+                m_min.y -= m_max.y - m_min.y;+            else if(point.y > m_max.y)+                m_max.y += m_max.y - m_min.y;+            break;+        case NodeSplitDirection::Y:+            m_rootDir = NodeSplitDirection::X;+            point.x < m_min.x ? m_nodes[newRoot].setChildren(newLeaf, m_root)+                              : m_nodes[newRoot].setChildren(m_root, newLeaf);+            if(point.x < m_min.x)+                m_min.x -= m_max.x - m_min.x;+            else if(point.x > m_max.x)+                m_max.x += m_max.x - m_min.x;+            break;+        }+        m_root = newRoot;+    }++    /// Calculate root's box enclosing all root points+    void initializeRootBox(const std::vector<point_type>& points)+    {+        const point_data_vec& data = m_nodes[m_root].data;+        m_min = points[data.front()];+        m_max = m_min;+        for(pd_cit it = data.begin(); it != data.end(); ++it)+        {+            const point_type& p = points[*it];+            m_min = point_type::make(+                std::min(m_min.x, p.x), std::min(m_min.y, p.y));+            m_max = point_type::make(+                std::max(m_max.x, p.x), std::max(m_max.y, p.y));+        }+        // Make sure bounding box does not have a zero size by adding padding:+        // zero-size bounding box cannot be extended properly+        const TCoordType padding(1);+        if(m_min.x == m_max.x)+        {+            m_min.x -= padding;+            m_max.x += padding;+        }+        if(m_min.y == m_max.y)+        {+            m_min.y -= padding;+            m_max.y += padding;+        }+        m_isRootBoxInitialized = true;+    }++private:+    node_index m_root;+    std::vector<Node> m_nodes;+    NodeSplitDirection::Enum m_rootDir;+    point_type m_min;+    point_type m_max;+    bool m_isRootBoxInitialized;++    // used for nearest query+    struct NearestTask+    {+        node_index node;+        point_type min, max;+        NodeSplitDirection::Enum dir;+        coord_type distSq;+        NearestTask()+        {}+        NearestTask(+            const node_index node,+            const point_type& min,+            const point_type& max,+            const NodeSplitDirection::Enum dir,+            const coord_type distSq)+            : node(node)+            , min(min)+            , max(max)+            , dir(dir)+            , distSq(distSq)+        {}+    };+    // allocated in class (not in the 'nearest' method) for better performance+    mutable std::vector<NearestTask> m_tasksStack;+};++} // namespace KDTree++#endif // header guard

@@ -0,0 +1,71 @@+/* 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/. */++/**+ * @file+ * Adapter between for KDTree and CDT+ */++#ifndef CDT_POINTKDTREE_H+#define CDT_POINTKDTREE_H++#include "CDTUtils.h"+#include "KDTree.h"++namespace CDT+{++/// KD-tree holding points+template <+    typename TCoordType,+    size_t NumVerticesInLeaf = 32,+    size_t InitialStackDepth = 32,+    size_t StackDepthIncrement = 32>+class LocatorKDTree+{+public:+    /// Initialize KD-tree with points+    void initialize(const std::vector<V2d<TCoordType> >& points)+    {+        typedef V2d<TCoordType> V2d_t;+        V2d_t min = points.front();+        V2d_t max = min;+        typedef typename std::vector<V2d_t>::const_iterator Cit;+        for(Cit it = points.begin(); it != points.end(); ++it)+        {+            min = V2d_t::make(std::min(min.x, it->x), std::min(min.y, it->y));+            max = V2d_t::make(std::max(max.x, it->x), std::max(max.y, it->y));+        }+        m_kdTree = KDTree_t(min, max);+        for(VertInd i = 0; i < points.size(); ++i)+        {+            m_kdTree.insert(i, points);+        }+    }+    /// Add point to KD-tree+    void addPoint(const VertInd i, const std::vector<V2d<TCoordType> >& points)+    {+        m_kdTree.insert(i, points);+    }+    /// Find nearest point using R-tree+    VertInd nearPoint(+        const V2d<TCoordType>& pos,+        const std::vector<V2d<TCoordType> >& points) const+    {+        return m_kdTree.nearest(pos, points).second;+    }++private:+    typedef KDTree::KDTree<+        TCoordType,+        NumVerticesInLeaf,+        InitialStackDepth,+        StackDepthIncrement>+        KDTree_t;+    KDTree_t m_kdTree;+};++} // namespace CDT++#endif

@@ -0,0 +1,617 @@+/* 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/. */++/**+ * @file+ * Triangulation class+ */++#ifndef CDT_vW1vZ0lO8rS4gY4uI4fB+#define CDT_vW1vZ0lO8rS4gY4uI4fB++#include "CDTUtils.h"+#include "LocatorKDTree.h"++#include <algorithm>+#include <cstdlib>+#include <iterator>+#include <stack>+#include <stdexcept>+#include <utility>+#include <vector>++/// Namespace containing triangulation functionality+namespace CDT+{++/// @addtogroup API+/// @{++/**+ * Enum of strategies specifying order in which a range of vertices is inserted+ * @note VertexInsertionOrder::Randomized will only randomize order of+ * inserting in triangulation, vertex indices will be preserved as they were+ * specified in the final triangulation+ */+struct CDT_EXPORT VertexInsertionOrder+{+    /**+     * The Enum itself+     * @note needed to pre c++11 compilers that don't support 'class enum'+     */+    enum Enum+    {+        Randomized, ///< vertices will be inserted in random order+        AsProvided, ///< vertices will be inserted in the same order as provided+    };+};++/// Enum of what type of geometry used to embed triangulation into+struct CDT_EXPORT SuperGeometryType+{+    /**+     * The Enum itself+     * @note needed to pre c++11 compilers that don't support 'class enum'+     */+    enum Enum+    {+        SuperTriangle, ///< conventional super-triangle+        Custom,        ///< user-specified custom geometry (e.g., grid)+    };+};++/**+ * Enum of strategies for treating intersecting constraint edges+ */+struct CDT_EXPORT IntersectingConstraintEdges+{+    /**+     * The Enum itself+     * @note needed to pre c++11 compilers that don't support 'class enum'+     */+    enum Enum+    {+        Ignore,  ///< constraint edge intersections are not checked+        Resolve, ///< constraint edge intersections are resolved+    };+};++/**+ * Type used for storing layer depths for triangles+ * @note LayerDepth should support 60K+ layers, which could be to much or+ * too little for some use cases. Feel free to re-define this typedef.+ */+typedef unsigned short LayerDepth;+typedef LayerDepth BoundaryOverlapCount;++/// Triangles by vertex index+typedef std::vector<TriIndVec> VerticesTriangles;++/**+ * @defgroup Triangulation Triangulation Class+ * Class performing triangulations.+ */+/// @{++/**+ * Data structure representing a 2D constrained Delaunay triangulation+ *+ * @tparam T type of vertex coordinates (e.g., float, double)+ * @tparam TNearPointLocator class providing locating near point for efficiently+ * inserting new points. Provides methods: 'addPoint(vPos, iV)' and+ * 'nearPoint(vPos) -> iV'+ */+template <typename T, typename TNearPointLocator = LocatorKDTree<T> >+class CDT_EXPORT Triangulation+{+public:+    typedef std::vector<V2d<T> > V2dVec; ///< Vertices vector+    V2dVec vertices;                     ///< triangulation's vertices+    TriangleVec triangles;               ///< triangulation's triangles+    EdgeUSet fixedEdges; ///< triangulation's constraints (fixed edges)+    /**+     * triangles adjacent to each vertex+     * @note will be reset to empty when super-triangle is removed and+     * triangulation is finalized. To re-calculate adjacent triangles use+     * CDT::calculateTrianglesByVertex helper+     */+    VerticesTriangles vertTris;++    /** Stores count of overlapping boundaries for a fixed edge. If no entry is+     * present for an edge: no boundaries overlap.+     * @note map only has entries for fixed for edges that represent overlapping+     * boundaries+     * @note needed for handling depth calculations and hole-removel in case of+     * overlapping boundaries+     */+    unordered_map<Edge, BoundaryOverlapCount> overlapCount;++    /** Stores list of original edges represented by a given fixed edge+     * @note map only has entries for edges where multiple original fixed edges+     * overlap or where a fixed edge is a part of original edge created by+     * conforming Delaunay triangulation vertex insertion+     */+    unordered_map<Edge, EdgeVec> pieceToOriginals;++    /*____ API _____*/+    /// Default constructor+    Triangulation();+    /**+     * Constructor+     * @param vertexInsertionOrder strategy used for ordering vertex insertions+     */+    Triangulation(VertexInsertionOrder::Enum vertexInsertionOrder);+    /**+     * Constructor+     * @param vertexInsertionOrder strategy used for ordering vertex insertions+     * @param intersectingEdgesStrategy strategy for treating intersecting+     * constraint edges+     * @param minDistToConstraintEdge distance within which point is considered+     * to be lying on a constraint edge. Used when adding constraints to the+     * triangulation.+     */+    Triangulation(+        VertexInsertionOrder::Enum vertexInsertionOrder,+        IntersectingConstraintEdges::Enum intersectingEdgesStrategy,+        T minDistToConstraintEdge);+    /**+     * Constructor+     * @param vertexInsertionOrder strategy used for ordering vertex insertions+     * @param nearPtLocator class providing locating near point for efficiently+     * inserting new points+     * @param intersectingEdgesStrategy strategy for treating intersecting+     * constraint edges+     * @param minDistToConstraintEdge distance within which point is considered+     * to be lying on a constraint edge. Used when adding constraints to the+     * triangulation.+     */+    Triangulation(+        VertexInsertionOrder::Enum vertexInsertionOrder,+        const TNearPointLocator& nearPtLocator,+        IntersectingConstraintEdges::Enum intersectingEdgesStrategy,+        T minDistToConstraintEdge);+    /**+     * Insert custom point-types specified by iterator range and X/Y-getters+     * @tparam TVertexIter iterator that dereferences to custom point type+     * @tparam TGetVertexCoordX function object getting x coordinate from+     * vertex. Getter signature: const TVertexIter::value_type& -> T+     * @tparam TGetVertexCoordY function object getting y coordinate from+     * vertex. Getter signature: const TVertexIter::value_type& -> T+     * @param first beginning of the range of vertices to add+     * @param last end of the range of vertices to add+     * @param getX getter of X-coordinate+     * @param getY getter of Y-coordinate+     */+    template <+        typename TVertexIter,+        typename TGetVertexCoordX,+        typename TGetVertexCoordY>+    void insertVertices(+        TVertexIter first,+        TVertexIter last,+        TGetVertexCoordX getX,+        TGetVertexCoordY getY);+    /**+     * Insert vertices into triangulation+     * @param vertices vector of vertices to insert+     */+    void insertVertices(const std::vector<V2d<T> >& vertices);+    /**+     * Insert constraints (custom-type fixed edges) into triangulation+     * @note Each fixed edge is inserted by deleting the triangles it crosses,+     * followed by the triangulation of the polygons on each side of the edge.+     * <b> No new vertices are inserted.</b>+     * @note If some edge appears more than once in the input this means that+     * multiple boundaries overlap at the edge and impacts how hole detection+     * algorithm of Triangulation::eraseOuterTrianglesAndHoles works.+     * <b>Make sure there are no erroneous duplicates.</b>+     * @tparam TEdgeIter iterator that dereferences to custom edge type+     * @tparam TGetEdgeVertexStart function object getting start vertex index+     * from an edge.+     * Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+     * @tparam TGetEdgeVertexEnd function object getting end vertex index from+     * an edge. Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+     * @param first beginning of the range of edges to add+     * @param last end of the range of edges to add+     * @param getStart getter of edge start vertex index+     * @param getEnd getter of edge end vertex index+     */+    template <+        typename TEdgeIter,+        typename TGetEdgeVertexStart,+        typename TGetEdgeVertexEnd>+    void insertEdges(+        TEdgeIter first,+        TEdgeIter last,+        TGetEdgeVertexStart getStart,+        TGetEdgeVertexEnd getEnd);+    /**+     * Insert constraint edges into triangulation+     * @note Each fixed edge is inserted by deleting the triangles it crosses,+     * followed by the triangulation of the polygons on each side of the edge.+     * <b> No new vertices are inserted.</b>+     * @note If some edge appears more than once in the input this means that+     * multiple boundaries overlap at the edge and impacts how hole detection+     * algorithm of Triangulation::eraseOuterTrianglesAndHoles works.+     * <b>Make sure there are no erroneous duplicates.</b>+     * @tparam edges constraint edges+     */+    void insertEdges(const std::vector<Edge>& edges);+    /**+     * Ensure that triangulation conforms to constraints (fixed edges)+     * @note For each fixed edge that is not present in the triangulation its+     * midpoint is recursively added until the original edge is represented by a+     * sequence of its pieces. <b> New vertices are inserted.</b>+     * @note If some edge appears more than once the input this+     * means that multiple boundaries overlap at the edge and impacts how hole+     * detection algorithm of Triangulation::eraseOuterTrianglesAndHoles works.+     * <b>Make sure there are no erroneous duplicates.</b>+     * @tparam TEdgeIter iterator that dereferences to custom edge type+     * @tparam TGetEdgeVertexStart function object getting start vertex index+     * from an edge.+     * Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+     * @tparam TGetEdgeVertexEnd function object getting end vertex index from+     * an edge. Getter signature: const TEdgeIter::value_type& -> CDT::VertInd+     * @param first beginning of the range of edges to add+     * @param last end of the range of edges to add+     * @param getStart getter of edge start vertex index+     * @param getEnd getter of edge end vertex index+     */+    template <+        typename TEdgeIter,+        typename TGetEdgeVertexStart,+        typename TGetEdgeVertexEnd>+    void conformToEdges(+        TEdgeIter first,+        TEdgeIter last,+        TGetEdgeVertexStart getStart,+        TGetEdgeVertexEnd getEnd);+    /**+     * Ensure that triangulation conforms to constraints (fixed edges)+     * @note For each fixed edge that is not present in the triangulation its+     * midpoint is recursively added until the original edge is represented by a+     * sequence of its pieces. <b> New vertices are inserted.</b>+     * @note If some edge appears more than once the input this+     * means that multiple boundaries overlap at the edge and impacts how hole+     * detection algorithm of Triangulation::eraseOuterTrianglesAndHoles works.+     * <b>Make sure there are no erroneous duplicates.</b>+     * @tparam edges edges to conform to+     */+    void conformToEdges(const std::vector<Edge>& edges);+    /**+     * Erase triangles adjacent to super triangle+     *+     * @note does nothing if custom geometry is used+     */+    void eraseSuperTriangle();+    /// Erase triangles outside of constrained boundary using growing+    void eraseOuterTriangles();+    /**+     * Erase triangles outside of constrained boundary and auto-detected holes+     *+     * @note detecting holes relies on layer peeling based on layer depth+     * @note supports overlapping or touching boundaries+     */+    void eraseOuterTrianglesAndHoles();+    /**+     * Call this method after directly setting custom super-geometry via+     * vertices and triangles members+     */+    void initializedWithCustomSuperGeometry();++    /**+     * Check if the triangulation was finalized with `erase...` method and+     * super-triangle was removed.+     * @return true if triangulation is finalized, false otherwise+     */+    bool isFinalized() const;++    /**+     * Calculate depth of each triangle in constraint triangulation. Supports+     * overlapping boundaries.+     *+     * Perform depth peeling from super triangle to outermost boundary,+     * then to next boundary and so on until all triangles are traversed.@n+     * For example depth is:+     *  - 0 for triangles outside outermost boundary+     *  - 1 for triangles inside boundary but outside hole+     *  - 2 for triangles in hole+     *  - 3 for triangles in island and so on...+     * @return vector where element at index i stores depth of i-th triangle+     */+    std::vector<LayerDepth> calculateTriangleDepths() const;++    /**+     * @defgroup Advanced Advanced Triangulation Methods+     * Advanced methods for manually modifying the triangulation from+     * outside. Please only use them when you know what you are doing.+     */+    /// @{++    /**+     * Flip an edge between two triangle.+     * @note Advanced method for manually modifying the triangulation from+     * outside. Please call it when you know what you are doing.+     * @param iT first triangle+     * @param iTopo second triangle++     */+    void flipEdge(TriInd iT, TriInd iTopo);++    /**+     * Remove triangles with specified indices.+     * Adjust internal triangulation state accordingly.+     * @param removedTriangles indices of triangles to remove+     */+    void removeTriangles(const TriIndUSet& removedTriangles);+    /// @}++private:+    /*____ Detail __*/+    void addSuperTriangle(const Box2d<T>& box);+    void addNewVertex(const V2d<T>& pos, const TriIndVec& tris);+    void insertVertex(VertInd iVert);+    void ensureDelaunayByEdgeFlips(+        const V2d<T>& v,+        VertInd iVert,+        std::stack<TriInd>& triStack);+    /// Flip fixed edges and return a list of flipped fixed edges+    std::vector<Edge> insertVertex_FlipFixedEdges(VertInd iVert);+    /**+     * Insert an edge into constraint Delaunay triangulation+     * @param edge edge to insert+     * @param originalEdge original edge inserted edge is part of+     */+    void insertEdge(Edge edge, Edge originalEdge);+    /**+     * Conform Delaunay triangulation to a fixed edge by recursively inserting+     * mid point of the edge and then conforming to its halves+     * @param edge fixed edge to conform to+     * @param originalEdges original edges that new edge is piece of+     * @param overlaps count of overlapping boundaries at the edge. Only used+     * when re-introducing edge with overlaps > 0+     * @param orientationTolerance tolerance for orient2d predicate,+     * values [-tolerance,+tolerance] are considered as 0.+     */+    void conformToEdge(+        Edge edge,+        EdgeVec originalEdges,+        BoundaryOverlapCount overlaps);+    tuple<TriInd, VertInd, VertInd> intersectedTriangle(+        VertInd iA,+        const std::vector<TriInd>& candidates,+        const V2d<T>& a,+        const V2d<T>& b,+        T orientationTolerance = T(0)) const;+    /// Returns indices of three resulting triangles+    std::stack<TriInd> insertPointInTriangle(VertInd v, TriInd iT);+    /// Returns indices of four resulting triangles+    std::stack<TriInd> insertPointOnEdge(VertInd v, TriInd iT1, TriInd iT2);+    array<TriInd, 2> trianglesAt(const V2d<T>& pos) const;+    array<TriInd, 2> walkingSearchTrianglesAt(const V2d<T>& pos) const;+    TriInd walkTriangles(VertInd startVertex, const V2d<T>& pos) const;+    bool isFlipNeeded(+        const V2d<T>& v,+        VertInd iV,+        VertInd iV1,+        VertInd iV2,+        VertInd iV3) const;+    bool+    isFlipNeeded(const V2d<T>& v, TriInd iT, TriInd iTopo, VertInd iVert) const;+    void changeNeighbor(TriInd iT, TriInd oldNeighbor, TriInd newNeighbor);+    void changeNeighbor(+        TriInd iT,+        VertInd iVedge1,+        VertInd iVedge2,+        TriInd newNeighbor);+    void addAdjacentTriangle(VertInd iVertex, TriInd iTriangle);+    void+    addAdjacentTriangles(VertInd iVertex, TriInd iT1, TriInd iT2, TriInd iT3);+    void addAdjacentTriangles(+        VertInd iVertex,+        TriInd iT1,+        TriInd iT2,+        TriInd iT3,+        TriInd iT4);+    void removeAdjacentTriangle(VertInd iVertex, TriInd iTriangle);+    TriInd triangulatePseudopolygon(+        VertInd ia,+        VertInd ib,+        std::vector<VertInd>::const_iterator pointsFirst,+        std::vector<VertInd>::const_iterator pointsLast);+    VertInd findDelaunayPoint(+        VertInd ia,+        VertInd ib,+        std::vector<VertInd>::const_iterator pointsFirst,+        std::vector<VertInd>::const_iterator pointsLast) const;+    TriInd pseudopolyOuterTriangle(VertInd ia, VertInd ib) const;+    TriInd addTriangle(const Triangle& t); // note: invalidates iterators!+    TriInd addTriangle(); // note: invalidates triangle iterators!+    /**+     * Remove super-triangle (if used) and triangles with specified indices.+     * Adjust internal triangulation state accordingly.+     * @removedTriangles indices of triangles to remove+     */+    void finalizeTriangulation(const TriIndUSet& removedTriangles);+    TriIndUSet growToBoundary(std::stack<TriInd> seeds) const;+    void fixEdge(const Edge& edge, BoundaryOverlapCount overlaps);+    void fixEdge(const Edge& edge);+    void fixEdge(const Edge& edge, const Edge& originalEdge);+    /**+     * Flag triangle as dummy+     * @note Advanced method for manually modifying the triangulation from+     * outside. Please call it when you know what you are doing.+     * @param iT index of a triangle to flag+     */+    void makeDummy(TriInd iT);+    /**+     * Erase all dummy triangles+     * @note Advanced method for manually modifying the triangulation from+     * outside. Please call it when you know what you are doing.+     */+    void eraseDummies();+    /**+     * Depth-peel a layer in triangulation, used when calculating triangle+     * depths+     *+     * It takes starting seed triangles, traverses neighboring triangles, and+     * assigns given layer depth to the traversed triangles. Traversal is+     * blocked by constraint edges. Triangles behind constraint edges are+     * recorded as seeds of next layer and returned from the function.+     *+     * @param seeds indices of seed triangles+     * @param layerDepth current layer's depth to mark triangles with+     * @param[in, out] triDepths depths of triangles+     * @return triangles of the deeper layers that are adjacent to the peeled+     * layer. To be used as seeds when peeling deeper layers.+     */+    unordered_map<TriInd, LayerDepth> peelLayer(+        std::stack<TriInd> seeds,+        LayerDepth layerDepth,+        std::vector<LayerDepth>& triDepths) const;++    std::vector<TriInd> m_dummyTris;+    TNearPointLocator m_nearPtLocator;+    std::size_t m_nTargetVerts;+    SuperGeometryType::Enum m_superGeomType;+    VertexInsertionOrder::Enum m_vertexInsertionOrder;+    IntersectingConstraintEdges::Enum m_intersectingEdgesStrategy;+    T m_minDistToConstraintEdge;+};++/// @}+/// @}++namespace detail+{++static mt19937 randGenerator(9001);++template <class RandomIt>+void random_shuffle(RandomIt first, RandomIt last)+{+    typename std::iterator_traits<RandomIt>::difference_type i, n;+    n = last - first;+    for(i = n - 1; i > 0; --i)+    {+        std::swap(first[i], first[randGenerator() % (i + 1)]);+    }+}++} // 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 std::runtime_error(+            "Triangulation was finalized with 'erase...' method. Inserting new "+            "vertices is not possible");+    }+    detail::randGenerator.seed(9001); // ensure deterministic behavior+    if(vertices.empty())+    {+        addSuperTriangle(envelopBox<T>(first, last, getX, getY));+    }++    const std::size_t nExistingVerts = vertices.size();++    vertices.reserve(nExistingVerts + std::distance(first, last));+    for(TVertexIter it = first; it != last; ++it)+        addNewVertex(V2d<T>::make(getX(*it), getY(*it)), TriIndVec());++    switch(m_vertexInsertionOrder)+    {+    case VertexInsertionOrder::AsProvided:+        for(TVertexIter it = first; it != last; ++it)+            insertVertex(VertInd(nExistingVerts + std::distance(first, it)));+        break;+    case VertexInsertionOrder::Randomized:+        std::vector<VertInd> ii(std::distance(first, last));+        typedef std::vector<VertInd>::iterator Iter;+        VertInd value = static_cast<VertInd>(nExistingVerts);+        for(Iter it = ii.begin(); it != ii.end(); ++it, ++value)+            *it = value;+        detail::random_shuffle(ii.begin(), ii.end());+        for(Iter it = ii.begin(); it != ii.end(); ++it)+            insertVertex(*it);+        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 std::runtime_error(+            "Triangulation was finalized with 'erase...' method. Inserting new "+            "edges is not possible");+    }+    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);+    }+    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 std::runtime_error(+            "Triangulation was finalized with 'erase...' method. Conforming to "+            "new edges is not possible");+    }+    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);+    }+    eraseDummies();+}++} // namespace CDT++#ifndef CDT_USE_AS_COMPILED_LIBRARY+#include "Triangulation.hpp"+#endif++#endif // header-guard

@@ -0,0 +1,1558 @@+/* 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/. */++/**+ * @file+ * Triangulation class - implementation+ */++#include "Triangulation.h"++#include <algorithm>+#include <cassert>+#include <deque>+#include <stdexcept>++namespace CDT+{++typedef std::deque<TriInd> TriDeque;++namespace detail+{++/// Needed for c++03 compatibility (no uniform initialization available)+template <typename T>+array<T, 3> arr3(const T& v0, const T& v1, const T& v2)+{+    const array<T, 3> out = {v0, v1, v2};+    return out;+}++namespace defaults+{++const std::size_t nTargetVerts = 0;+const SuperGeometryType::Enum superGeomType = SuperGeometryType::SuperTriangle;+const VertexInsertionOrder::Enum vertexInsertionOrder =+    VertexInsertionOrder::Randomized;+const IntersectingConstraintEdges::Enum intersectingEdgesStrategy =+    IntersectingConstraintEdges::Ignore;+const float minDistToConstraintEdge(0);++} // namespace defaults++} // namespace detail++template <typename T, typename TNearPointLocator>+Triangulation<T, TNearPointLocator>::Triangulation()+    : m_nTargetVerts(detail::defaults::nTargetVerts)+    , m_superGeomType(detail::defaults::superGeomType)+    , m_vertexInsertionOrder(detail::defaults::vertexInsertionOrder)+    , m_intersectingEdgesStrategy(detail::defaults::intersectingEdgesStrategy)+    , m_minDistToConstraintEdge(detail::defaults::minDistToConstraintEdge)+{}++template <typename T, typename TNearPointLocator>+Triangulation<T, TNearPointLocator>::Triangulation(+    const VertexInsertionOrder::Enum vertexInsertionOrder)+    : m_nTargetVerts(detail::defaults::nTargetVerts)+    , m_superGeomType(detail::defaults::superGeomType)+    , m_vertexInsertionOrder(vertexInsertionOrder)+    , m_intersectingEdgesStrategy(detail::defaults::intersectingEdgesStrategy)+    , m_minDistToConstraintEdge(detail::defaults::minDistToConstraintEdge)+{}++template <typename T, typename TNearPointLocator>+Triangulation<T, TNearPointLocator>::Triangulation(+    const VertexInsertionOrder::Enum vertexInsertionOrder,+    const IntersectingConstraintEdges::Enum intersectingEdgesStrategy,+    const T minDistToConstraintEdge)+    : m_nTargetVerts(detail::defaults::nTargetVerts)+    , m_superGeomType(detail::defaults::superGeomType)+    , m_vertexInsertionOrder(vertexInsertionOrder)+    , m_intersectingEdgesStrategy(intersectingEdgesStrategy)+    , m_minDistToConstraintEdge(minDistToConstraintEdge)+{}++template <typename T, typename TNearPointLocator>+Triangulation<T, TNearPointLocator>::Triangulation(+    const VertexInsertionOrder::Enum vertexInsertionOrder,+    const TNearPointLocator& nearPtLocator,+    const IntersectingConstraintEdges::Enum intersectingEdgesStrategy,+    const T minDistToConstraintEdge)+    : m_nTargetVerts(detail::defaults::nTargetVerts)+    , m_nearPtLocator(nearPtLocator)+    , m_superGeomType(detail::defaults::superGeomType)+    , m_vertexInsertionOrder(vertexInsertionOrder)+    , m_intersectingEdgesStrategy(intersectingEdgesStrategy)+    , m_minDistToConstraintEdge(minDistToConstraintEdge)+{}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::changeNeighbor(+    const TriInd iT,+    const VertInd iVedge1,+    const VertInd iVedge2,+    const TriInd newNeighbor)+{+    Triangle& t = triangles[iT];+    t.neighbors[opposedTriangleInd(t, iVedge1, iVedge2)] = newNeighbor;+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::eraseDummies()+{+    if(m_dummyTris.empty())+        return;+    const TriIndUSet dummySet(m_dummyTris.begin(), m_dummyTris.end());+    TriIndUMap triIndMap;+    triIndMap[noNeighbor] = noNeighbor;+    for(TriInd iT(0), iTnew(0); iT < TriInd(triangles.size()); ++iT)+    {+        if(dummySet.count(iT))+            continue;+        triIndMap[iT] = iTnew;+        triangles[iTnew] = triangles[iT];+        iTnew++;+    }+    triangles.erase(triangles.end() - dummySet.size(), triangles.end());++    // remap adjacent triangle indices for vertices+    typedef typename VerticesTriangles::iterator VertTrisIt;+    for(VertTrisIt vTris = vertTris.begin(); vTris != vertTris.end(); ++vTris)+    {+        for(TriIndVec::iterator iT = vTris->begin(); iT != vTris->end(); ++iT)+            *iT = triIndMap[*iT];+    }+    // remap neighbor indices for triangles+    for(TriangleVec::iterator t = triangles.begin(); t != triangles.end(); ++t)+    {+        NeighborsArr3& nn = t->neighbors;+        for(NeighborsArr3::iterator iN = nn.begin(); iN != nn.end(); ++iN)+            *iN = triIndMap[*iN];+    }+    // clear dummy triangles+    m_dummyTris = std::vector<TriInd>();+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::eraseSuperTriangle()+{+    if(m_superGeomType != SuperGeometryType::SuperTriangle)+        return;+    // find triangles adjacent to super-triangle's vertices+    TriIndUSet toErase;+    toErase.reserve(+        vertTris[0].size() + vertTris[1].size() + vertTris[2].size());+    for(TriInd iT(0); iT < TriInd(triangles.size()); ++iT)+    {+        Triangle& t = triangles[iT];+        if(t.vertices[0] < 3 || t.vertices[1] < 3 || t.vertices[2] < 3)+            toErase.insert(iT);+    }+    finalizeTriangulation(toErase);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::eraseOuterTriangles()+{+    // make dummy triangles adjacent to super-triangle's vertices+    const std::stack<TriInd> seed(std::deque<TriInd>(1, vertTris[0].front()));+    const TriIndUSet toErase = growToBoundary(seed);+    finalizeTriangulation(toErase);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::eraseOuterTrianglesAndHoles()+{+    const std::vector<LayerDepth> triDepths = calculateTriangleDepths();+    TriIndUSet toErase;+    toErase.reserve(triangles.size());+    for(std::size_t iT = 0; iT != triangles.size(); ++iT)+    {+        if(triDepths[iT] % 2 == 0)+            toErase.insert(static_cast<TriInd>(iT));+    }+    finalizeTriangulation(toErase);+}++/// Remap removing super-triangle: subtract 3 from vertices+inline Edge RemapNoSuperTriangle(const Edge& e)+{+    return Edge(e.v1() - 3, e.v2() - 3);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::removeTriangles(+    const TriIndUSet& removedTriangles)+{+    if(removedTriangles.empty())+        return;+    // remove triangles and calculate triangle index mapping+    TriIndUMap triIndMap;+    for(TriInd iT(0), iTnew(0); iT < TriInd(triangles.size()); ++iT)+    {+        if(removedTriangles.count(iT))+            continue;+        triIndMap[iT] = iTnew;+        triangles[iTnew] = triangles[iT];+        iTnew++;+    }+    triangles.erase(triangles.end() - removedTriangles.size(), triangles.end());+    // adjust triangles' neighbors+    vertTris = VerticesTriangles();+    for(TriInd iT = 0; iT < triangles.size(); ++iT)+    {+        Triangle& t = triangles[iT];+        // update neighbors to account for removed triangles+        NeighborsArr3& nn = t.neighbors;+        for(NeighborsArr3::iterator n = nn.begin(); n != nn.end(); ++n)+        {+            if(removedTriangles.count(*n))+            {+                *n = noNeighbor;+            }+            else if(*n != noNeighbor)+            {+                *n = triIndMap[*n];+            }+        }+    }+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::finalizeTriangulation(+    const TriIndUSet& removedTriangles)+{+    eraseDummies();+    // remove super-triangle+    if(m_superGeomType == SuperGeometryType::SuperTriangle)+    {+        vertices.erase(vertices.begin(), vertices.begin() + 3);+        if(removedTriangles.empty())+            vertTris.erase(vertTris.begin(), vertTris.begin() + 3);+        // Edge re-mapping+        { // fixed edges+            EdgeUSet updatedFixedEdges;+            typedef CDT::EdgeUSet::const_iterator It;+            for(It e = fixedEdges.begin(); e != fixedEdges.end(); ++e)+            {+                updatedFixedEdges.insert(RemapNoSuperTriangle(*e));+            }+            fixedEdges = updatedFixedEdges;+        }+        { // overlap count+            unordered_map<Edge, BoundaryOverlapCount> updatedOverlapCount;+            typedef unordered_map<Edge, BoundaryOverlapCount>::const_iterator+                It;+            for(It it = overlapCount.begin(); it != overlapCount.end(); ++it)+            {+                updatedOverlapCount.insert(std::make_pair(+                    RemapNoSuperTriangle(it->first), it->second));+            }+            overlapCount = updatedOverlapCount;+        }+        { // split edges mapping+            unordered_map<Edge, EdgeVec> updatedPieceToOriginals;+            typedef unordered_map<Edge, EdgeVec>::const_iterator It;+            for(It it = pieceToOriginals.begin(); it != pieceToOriginals.end();+                ++it)+            {+                EdgeVec ee = it->second;+                for(EdgeVec::iterator eeIt = ee.begin(); eeIt != ee.end();+                    ++eeIt)+                {+                    *eeIt = RemapNoSuperTriangle(*eeIt);+                }+                updatedPieceToOriginals.insert(+                    std::make_pair(RemapNoSuperTriangle(it->first), ee));+            }+            pieceToOriginals = updatedPieceToOriginals;+        }+    }+    // remove other triangles+    removeTriangles(removedTriangles);+    // adjust triangle vertices: account for removed super-triangle+    if(m_superGeomType == SuperGeometryType::SuperTriangle)+    {+        for(TriangleVec::iterator t = triangles.begin(); t != triangles.end();+            ++t)+        {+            VerticesArr3& vv = t->vertices;+            for(VerticesArr3::iterator v = vv.begin(); v != vv.end(); ++v)+            {+                *v -= 3;+            }+        }+    }+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::initializedWithCustomSuperGeometry()+{+    m_nearPtLocator.initialize(vertices);+    m_nTargetVerts = vertices.size();+    m_superGeomType = SuperGeometryType::Custom;+}++template <typename T, typename TNearPointLocator>+TriIndUSet Triangulation<T, TNearPointLocator>::growToBoundary(+    std::stack<TriInd> seeds) const+{+    TriIndUSet traversed;+    while(!seeds.empty())+    {+        const TriInd iT = seeds.top();+        seeds.pop();+        traversed.insert(iT);+        const Triangle& t = triangles[iT];+        for(Index i(0); i < Index(3); ++i)+        {+            const Edge opEdge(t.vertices[ccw(i)], t.vertices[cw(i)]);+            if(fixedEdges.count(opEdge))+                continue;+            const TriInd iN = t.neighbors[opoNbr(i)];+            if(iN != noNeighbor && traversed.count(iN) == 0)+                seeds.push(iN);+        }+    }+    return traversed;+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::makeDummy(const TriInd iT)+{+    const Triangle& t = triangles[iT];++    typedef VerticesArr3::const_iterator VCit;+    for(VCit iV = t.vertices.begin(); iV != t.vertices.end(); ++iV)+        removeAdjacentTriangle(*iV, iT);++    typedef NeighborsArr3::const_iterator NCit;+    for(NCit iTn = t.neighbors.begin(); iTn != t.neighbors.end(); ++iTn)+        changeNeighbor(*iTn, iT, noNeighbor);++    m_dummyTris.push_back(iT);+}++template <typename T, typename TNearPointLocator>+TriInd Triangulation<T, TNearPointLocator>::addTriangle(const Triangle& t)+{+    if(m_dummyTris.empty())+    {+        triangles.push_back(t);+        return TriInd(triangles.size() - 1);+    }+    const TriInd nxtDummy = m_dummyTris.back();+    m_dummyTris.pop_back();+    triangles[nxtDummy] = t;+    return nxtDummy;+}++template <typename T, typename TNearPointLocator>+TriInd Triangulation<T, TNearPointLocator>::addTriangle()+{+    if(m_dummyTris.empty())+    {+        const Triangle dummy = {+            {noVertex, noVertex, noVertex},+            {noNeighbor, noNeighbor, noNeighbor}};+        triangles.push_back(dummy);+        return TriInd(triangles.size() - 1);+    }+    const TriInd nxtDummy = m_dummyTris.back();+    m_dummyTris.pop_back();+    return nxtDummy;+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::insertEdges(+    const std::vector<Edge>& edges)+{+    insertEdges(edges.begin(), edges.end(), edge_get_v1, edge_get_v2);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::conformToEdges(+    const std::vector<Edge>& edges)+{+    conformToEdges(edges.begin(), edges.end(), edge_get_v1, edge_get_v2);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::fixEdge(const Edge& edge)+{+    if(!fixedEdges.insert(edge).second)+    {+        ++overlapCount[edge]; // if edge is already fixed increment the counter+    }+}++namespace detail+{++// add element to 'to' if not already in 'to'+template <typename T, typename Allocator1>+void insert_unique(std::vector<T, Allocator1>& to, const T& elem)+{+    if(std::find(to.begin(), to.end(), elem) == to.end())+    {+        to.push_back(elem);+    }+}++// add elements of 'from' that are not present in 'to' to 'to'+template <typename T, typename Allocator1, typename Allocator2>+void insert_unique(+    std::vector<T, Allocator1>& to,+    const std::vector<T, Allocator2>& from)+{+    typedef typename std::vector<T, Allocator2>::const_iterator Cit;+    to.reserve(to.size() + from.size());+    for(Cit cit = from.begin(); cit != from.end(); ++cit)+    {+        insert_unique(to, *cit);+    }+}++} // namespace detail++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::fixEdge(+    const Edge& edge,+    const Edge& originalEdge)+{+    fixEdge(edge);+    if(edge != originalEdge)+        detail::insert_unique(pieceToOriginals[edge], originalEdge);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::fixEdge(+    const Edge& edge,+    const BoundaryOverlapCount overlaps)+{+    fixedEdges.insert(edge);+    overlapCount[edge] = overlaps; // override overlap counter+}++namespace detail+{++template <typename T>+T lerp(const T& a, const T& b, const T t)+{+    return (T(1) - t) * a + t * b;+}++// Precondition: ab and cd intersect normally+template <typename T>+V2d<T> intersectionPosition(+    const V2d<T>& a,+    const V2d<T>& b,+    const V2d<T>& c,+    const V2d<T>& d)+{+    using namespace predicates::adaptive;+    // interpolate point on the shorter segment+    if(distanceSquared(a, b) < distanceSquared(c, d))+    {+        const T a_cd = orient2d(c.x, c.y, d.x, d.y, a.x, a.y);+        const T b_cd = orient2d(c.x, c.y, d.x, d.y, b.x, b.y);+        const T t = a_cd / (a_cd - b_cd);+        return V2d<T>::make(lerp(a.x, b.x, t), lerp(a.y, b.y, t));+    }+    else+    {+        const T c_ab = orient2d(a.x, a.y, b.x, b.y, c.x, c.y);+        const T d_ab = orient2d(a.x, a.y, b.x, b.y, d.x, d.y);+        const T t = c_ab / (c_ab - d_ab);+        return V2d<T>::make(lerp(c.x, d.x, t), lerp(c.y, d.y, t));+    }+}++} // namespace detail++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::insertEdge(+    const Edge edge,+    const Edge originalEdge)+{+    const VertInd iA = edge.v1();+    VertInd iB = edge.v2();+    if(iA == iB) // edge connects a vertex to itself+        return;+    const TriIndVec& aTris = vertTris[iA];+    const TriIndVec& bTris = vertTris[iB];+    const V2d<T>& a = vertices[iA];+    const V2d<T>& b = vertices[iB];+    if(verticesShareEdge(aTris, bTris))+    {+        fixEdge(edge, originalEdge);+        return;+    }++    const T distanceTolerance =+        m_minDistToConstraintEdge == T(0)+            ? T(0)+            : m_minDistToConstraintEdge * distance(a, b);++    TriInd iT;+    VertInd iVleft, iVright;+    tie(iT, iVleft, iVright) =+        intersectedTriangle(iA, aTris, a, b, distanceTolerance);+    // if one of the triangle vertices is on the edge, move edge start+    if(iT == noNeighbor)+    {+        const Edge edgePart(iA, iVleft);+        fixEdge(edgePart, originalEdge);+        return insertEdge(Edge(iVleft, iB), originalEdge);+    }+    std::vector<TriInd> intersected(1, iT);+    std::vector<VertInd> ptsLeft(1, iVleft);+    std::vector<VertInd> ptsRight(1, iVright);+    VertInd iV = iA;+    Triangle t = triangles[iT];+    while(std::find(t.vertices.begin(), t.vertices.end(), iB) ==+          t.vertices.end())+    {+        const TriInd iTopo = opposedTriangle(t, iV);+        const Triangle& tOpo = triangles[iTopo];+        const VertInd iVopo = opposedVertex(tOpo, iT);+        const V2d<T> vOpo = vertices[iVopo];++        // Resolve intersection between two constraint edges if needed+        if(m_intersectingEdgesStrategy ==+               IntersectingConstraintEdges::Resolve &&+           fixedEdges.count(Edge(iVleft, iVright)))+        {+            const VertInd iNewVert = static_cast<VertInd>(vertices.size());++            // split constraint edge that already exists in triangulation+            const Edge splitEdge(iVleft, iVright);+            const Edge half1(iVleft, iNewVert);+            const Edge half2(iNewVert, iVright);+            const BoundaryOverlapCount overlaps = overlapCount[splitEdge];+            // remove the edge that will be split+            fixedEdges.erase(splitEdge);+            overlapCount.erase(splitEdge);+            // add split edge's halves+            fixEdge(half1, overlaps);+            fixEdge(half2, overlaps);+            // maintain piece-to-original mapping+            EdgeVec newOriginals(1, splitEdge);+            const unordered_map<Edge, EdgeVec>::const_iterator originalsIt =+                pieceToOriginals.find(splitEdge);+            if(originalsIt != pieceToOriginals.end())+            { // edge being split was split before: pass-through originals+                newOriginals = originalsIt->second;+                pieceToOriginals.erase(originalsIt);+            }+            detail::insert_unique(pieceToOriginals[half1], newOriginals);+            detail::insert_unique(pieceToOriginals[half2], newOriginals);++            // add a new point at the intersection of two constraint edges+            const V2d<T> newV = detail::intersectionPosition(+                vertices[iA],+                vertices[iB],+                vertices[iVleft],+                vertices[iVright]);+            addNewVertex(newV, TriIndVec());+            std::stack<TriInd> triStack =+                insertPointOnEdge(iNewVert, iT, iTopo);+            ensureDelaunayByEdgeFlips(newV, iNewVert, triStack);+            // TODO: is it's possible to re-use pseudo-polygons+            //  for inserting [iA, iNewVert] edge half?+            insertEdge(Edge(iA, iNewVert), originalEdge);+            insertEdge(Edge(iNewVert, iB), originalEdge);+            return;+        }++        intersected.push_back(iTopo);+        iT = iTopo;+        t = triangles[iT];++        const PtLineLocation::Enum loc =+            locatePointLine(vOpo, a, b, distanceTolerance);+        if(loc == PtLineLocation::Left)+        {+            ptsLeft.push_back(iVopo);+            iV = iVleft;+            iVleft = iVopo;+        }+        else if(loc == PtLineLocation::Right)+        {+            ptsRight.push_back(iVopo);+            iV = iVright;+            iVright = iVopo;+        }+        else // encountered point on the edge+            iB = iVopo;+    }+    // Remove intersected triangles+    typedef std::vector<TriInd>::const_iterator TriIndCit;+    for(TriIndCit it = intersected.begin(); it != intersected.end(); ++it)+        makeDummy(*it);+    // Triangulate pseudo-polygons on both sides+    const TriInd iTleft =+        triangulatePseudopolygon(iA, iB, ptsLeft.begin(), ptsLeft.end());+    std::reverse(ptsRight.begin(), ptsRight.end());+    const TriInd iTright =+        triangulatePseudopolygon(iB, iA, ptsRight.begin(), ptsRight.end());+    changeNeighbor(iTleft, noNeighbor, iTright);+    changeNeighbor(iTright, noNeighbor, iTleft);++    if(iB != edge.v2()) // encountered point on the edge+    {+        // fix edge part+        const Edge edgePart(iA, iB);+        fixEdge(edgePart, originalEdge);+        return insertEdge(Edge(iB, edge.v2()), originalEdge);+    }+    else+    {+        fixEdge(edge, originalEdge);+    }+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::conformToEdge(+    const Edge edge,+    EdgeVec originalEdges,+    const BoundaryOverlapCount overlaps)+{+    const VertInd iA = edge.v1();+    VertInd iB = edge.v2();+    if(iA == iB) // edge connects a vertex to itself+        return;+    const TriIndVec& aTris = vertTris[iA];+    const TriIndVec& bTris = vertTris[iB];+    const V2d<T>& a = vertices[iA];+    const V2d<T>& b = vertices[iB];+    if(verticesShareEdge(aTris, bTris))+    {+        overlaps > 0 ? fixEdge(edge, overlaps) : fixEdge(edge);+        // avoid marking edge as a part of itself+        if(!originalEdges.empty() && edge != originalEdges.front())+        {+            detail::insert_unique(pieceToOriginals[edge], originalEdges);+        }+        return;+    }++    const T distanceTolerance =+        m_minDistToConstraintEdge == T(0)+            ? T(0)+            : m_minDistToConstraintEdge * distance(a, b);+    TriInd iT;+    VertInd iVleft, iVright;+    tie(iT, iVleft, iVright) =+        intersectedTriangle(iA, aTris, a, b, distanceTolerance);+    // if one of the triangle vertices is on the edge, move edge start+    if(iT == noNeighbor)+    {+        const Edge edgePart(iA, iVleft);+        overlaps > 0 ? fixEdge(edgePart, overlaps) : fixEdge(edgePart);+        detail::insert_unique(pieceToOriginals[edgePart], originalEdges);+        return conformToEdge(Edge(iVleft, iB), originalEdges, overlaps);+    }++    VertInd iV = iA;+    Triangle t = triangles[iT];+    while(std::find(t.vertices.begin(), t.vertices.end(), iB) ==+          t.vertices.end())+    {+        const TriInd iTopo = opposedTriangle(t, iV);+        const Triangle& tOpo = triangles[iTopo];+        const VertInd iVopo = opposedVertex(tOpo, iT);+        const V2d<T> vOpo = vertices[iVopo];++        // Resolve intersection between two constraint edges if needed+        if(m_intersectingEdgesStrategy ==+               IntersectingConstraintEdges::Resolve &&+           fixedEdges.count(Edge(iVleft, iVright)))+        {+            const VertInd iNewVert = static_cast<VertInd>(vertices.size());++            // split constraint edge that already exists in triangulation+            const Edge splitEdge(iVleft, iVright);+            const Edge half1(iVleft, iNewVert);+            const Edge half2(iNewVert, iVright);+            const BoundaryOverlapCount overlaps = overlapCount[splitEdge];+            // remove the edge that will be split+            fixedEdges.erase(splitEdge);+            overlapCount.erase(splitEdge);+            // add split edge's halves+            fixEdge(half1, overlaps);+            fixEdge(half2, overlaps);+            // maintain piece-to-original mapping+            EdgeVec newOriginals(1, splitEdge);+            const unordered_map<Edge, EdgeVec>::const_iterator originalsIt =+                pieceToOriginals.find(splitEdge);+            if(originalsIt != pieceToOriginals.end())+            { // edge being split was split before: pass-through originals+                newOriginals = originalsIt->second;+                pieceToOriginals.erase(originalsIt);+            }+            detail::insert_unique(pieceToOriginals[half1], newOriginals);+            detail::insert_unique(pieceToOriginals[half2], newOriginals);++            // add a new point at the intersection of two constraint edges+            const V2d<T> newV = detail::intersectionPosition(+                vertices[iA],+                vertices[iB],+                vertices[iVleft],+                vertices[iVright]);+            addNewVertex(newV, TriIndVec());+            std::stack<TriInd> triStack =+                insertPointOnEdge(iNewVert, iT, iTopo);+            ensureDelaunayByEdgeFlips(newV, iNewVert, triStack);+            conformToEdge(Edge(iA, iNewVert), originalEdges, overlaps);+            conformToEdge(Edge(iNewVert, iB), originalEdges, overlaps);+            return;+        }++        iT = iTopo;+        t = triangles[iT];++        const PtLineLocation::Enum loc =+            locatePointLine(vOpo, a, b, distanceTolerance);+        if(loc == PtLineLocation::Left)+        {+            iV = iVleft;+            iVleft = iVopo;+        }+        else if(loc == PtLineLocation::Right)+        {+            iV = iVright;+            iVright = iVopo;+        }+        else // encountered point on the edge+            iB = iVopo;+    }+    /**/++    // add mid-point to triangulation+    const VertInd iMid = static_cast<VertInd>(vertices.size());+    const V2d<T>& start = vertices[iA];+    const V2d<T>& end = vertices[iB];+    addNewVertex(+        V2d<T>::make((start.x + end.x) / T(2), (start.y + end.y) / T(2)),+        TriIndVec());+    const std::vector<Edge> flippedFixedEdges =+        insertVertex_FlipFixedEdges(iMid);++    conformToEdge(Edge(iA, iMid), originalEdges, overlaps);+    conformToEdge(Edge(iMid, iB), originalEdges, overlaps);+    // re-introduce fixed edges that were flipped+    // and make sure overlap count is preserved+    for(std::vector<Edge>::const_iterator it = flippedFixedEdges.begin();+        it != flippedFixedEdges.end();+        ++it)+    {+        fixedEdges.erase(*it);++        BoundaryOverlapCount prevOverlaps = 0;+        const unordered_map<Edge, BoundaryOverlapCount>::const_iterator+            overlapsIt = overlapCount.find(*it);+        if(overlapsIt != overlapCount.end())+        {+            prevOverlaps = overlapsIt->second;+            overlapCount.erase(overlapsIt);+        }+        // override overlapping boundaries count when re-inserting an edge+        EdgeVec prevOriginals(1, *it);+        const unordered_map<Edge, EdgeVec>::const_iterator originalsIt =+            pieceToOriginals.find(*it);+        if(originalsIt != pieceToOriginals.end())+        {+            prevOriginals = originalsIt->second;+        }+        conformToEdge(*it, prevOriginals, prevOverlaps);+    }+    if(iB != edge.v2())+        conformToEdge(Edge(iB, edge.v2()), originalEdges, overlaps);+}++/*!+ * Returns:+ *  - intersected triangle index+ *  - index of point on the left of the line+ *  - index of point on the right of the line+ * If left point is right on the line: no triangle is intersected:+ *  - triangle index is no-neighbor (invalid)+ *  - index of point on the line+ *  - index of point on the right of the line+ */+template <typename T, typename TNearPointLocator>+tuple<TriInd, VertInd, VertInd>+Triangulation<T, TNearPointLocator>::intersectedTriangle(+    const VertInd iA,+    const std::vector<TriInd>& candidates,+    const V2d<T>& a,+    const V2d<T>& b,+    const T orientationTolerance) const+{+    typedef std::vector<TriInd>::const_iterator TriIndCit;+    for(TriIndCit it = candidates.begin(); it != candidates.end(); ++it)+    {+        const TriInd iT = *it;+        const Triangle t = triangles[iT];+        const Index i = vertexInd(t, iA);+        const VertInd iP2 = t.vertices[ccw(i)];+        const T orientP2 = orient2D(vertices[iP2], a, b);+        const PtLineLocation::Enum locP2 = classifyOrientation(orientP2);+        if(locP2 == PtLineLocation::Right)+        {+            const VertInd iP1 = t.vertices[cw(i)];+            const T orientP1 = orient2D(vertices[iP1], a, b);+            const PtLineLocation::Enum locP1 = classifyOrientation(orientP1);+            if(locP1 == PtLineLocation::OnLine)+            {+                return make_tuple(noNeighbor, iP1, iP1);+            }+            if(locP1 == PtLineLocation::Left)+            {+                if(orientationTolerance)+                {+                    T closestOrient;+                    VertInd iClosestP;+                    if(std::abs(orientP1) <= std::abs(orientP2))+                    {+                        closestOrient = orientP1;+                        iClosestP = iP1;+                    }+                    else+                    {+                        closestOrient = orientP2;+                        iClosestP = iP2;+                    }+                    if(classifyOrientation(+                           closestOrient, orientationTolerance) ==+                       PtLineLocation::OnLine)+                    {+                        return make_tuple(noNeighbor, iClosestP, iClosestP);+                    }+                }+                return make_tuple(iT, iP1, iP2);+            }+        }+    }+    throw std::runtime_error("Could not find vertex triangle intersected by "+                             "edge. Note: can be caused by duplicate points.");+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::addSuperTriangle(const Box2d<T>& box)+{+    m_nTargetVerts = 3;+    m_superGeomType = SuperGeometryType::SuperTriangle;++    const V2d<T> center = {+        (box.min.x + box.max.x) / T(2), (box.min.y + box.max.y) / T(2)};+    const T w = box.max.x - box.min.x;+    const T h = box.max.y - box.min.y;+    T r = std::sqrt(w * w + h * h) / T(2); // incircle radius+    r *= T(1.1);+    const T R = T(2) * r;                        // excircle radius+    const T shiftX = R * std::sqrt(T(3)) / T(2); // R * cos(30 deg)+    const V2d<T> posV1 = {center.x - shiftX, center.y - r};+    const V2d<T> posV2 = {center.x + shiftX, center.y - r};+    const V2d<T> posV3 = {center.x, center.y + R};+    addNewVertex(posV1, TriIndVec(1, TriInd(0)));+    addNewVertex(posV2, TriIndVec(1, TriInd(0)));+    addNewVertex(posV3, TriIndVec(1, TriInd(0)));+    const Triangle superTri = {+        {VertInd(0), VertInd(1), VertInd(2)},+        {noNeighbor, noNeighbor, noNeighbor}};+    addTriangle(superTri);+    m_nearPtLocator.initialize(vertices);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::addNewVertex(+    const V2d<T>& pos,+    const TriIndVec& tris)+{+    vertices.push_back(pos);+    vertTris.push_back(tris);+}++template <typename T, typename TNearPointLocator>+std::vector<Edge>+Triangulation<T, TNearPointLocator>::insertVertex_FlipFixedEdges(+    const VertInd iVert)+{+    std::vector<Edge> flippedFixedEdges;++    const V2d<T>& v = vertices[iVert];+    array<TriInd, 2> trisAt = walkingSearchTrianglesAt(v);+    std::stack<TriInd> triStack =+        trisAt[1] == noNeighbor+            ? insertPointInTriangle(iVert, trisAt[0])+            : insertPointOnEdge(iVert, trisAt[0], trisAt[1]);+    while(!triStack.empty())+    {+        const TriInd iT = triStack.top();+        triStack.pop();++        const Triangle& t = triangles[iT];+        const TriInd iTopo = opposedTriangle(t, iVert);+        if(iTopo == noNeighbor)+            continue;++        /*+         *                       v3         original edge: (v1, v3)+         *                      /|\   flip-candidate edge: (v,  v2)+         *                    /  |  \+         *                  /    |    \+         *                /      |      \+         * new vertex--> v       |       v2+         *                \      |      /+         *                  \    |    /+         *                    \  |  /+         *                      \|/+         *                       v1+         */+        const Triangle& tOpo = triangles[iTopo];+        const Index i = opposedVertexInd(tOpo, iT);+        const VertInd iV2 = tOpo.vertices[i];+        const VertInd iV1 = tOpo.vertices[cw(i)];+        const VertInd iV3 = tOpo.vertices[ccw(i)];++        if(isFlipNeeded(v, iVert, iV1, iV2, iV3))+        {+            // if flipped edge is fixed, remember it+            const Edge flippedEdge(iV1, iV3);+            if(fixedEdges.count(flippedEdge))+                flippedFixedEdges.push_back(flippedEdge);++            flipEdge(iT, iTopo);+            triStack.push(iT);+            triStack.push(iTopo);+        }+    }++    m_nearPtLocator.addPoint(iVert, vertices);+    return flippedFixedEdges;+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::insertVertex(const VertInd iVert)+{+    const V2d<T>& v = vertices[iVert];+    array<TriInd, 2> trisAt = walkingSearchTrianglesAt(v);+    std::stack<TriInd> triStack =+        trisAt[1] == noNeighbor+            ? insertPointInTriangle(iVert, trisAt[0])+            : insertPointOnEdge(iVert, trisAt[0], trisAt[1]);+    ensureDelaunayByEdgeFlips(v, iVert, triStack);+    m_nearPtLocator.addPoint(iVert, vertices);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::ensureDelaunayByEdgeFlips(+    const V2d<T>& v,+    const VertInd iVert,+    std::stack<TriInd>& triStack)+{+    while(!triStack.empty())+    {+        const TriInd iT = triStack.top();+        triStack.pop();++        const Triangle& t = triangles[iT];+        const TriInd iTopo = opposedTriangle(t, iVert);+        if(iTopo == noNeighbor)+            continue;+        if(isFlipNeeded(v, iT, iTopo, iVert))+        {+            flipEdge(iT, iTopo);+            triStack.push(iT);+            triStack.push(iTopo);+        }+    }+}++/*!+ * Handles super-triangle vertices.+ * Super-tri points are not infinitely far and influence the input points+ * Three cases are possible:+ *  1.  If one of the opposed vertices is super-tri: no flip needed+ *  2.  One of the shared vertices is super-tri:+ *      check if on point is same side of line formed by non-super-tri+ * vertices as the non-super-tri shared vertex+ *  3.  None of the vertices are super-tri: normal circumcircle test+ */+/*+ *                       v3         original edge: (v1, v3)+ *                      /|\   flip-candidate edge: (v,  v2)+ *                    /  |  \+ *                  /    |    \+ *                /      |      \+ * new vertex--> v       |       v2+ *                \      |      /+ *                  \    |    /+ *                    \  |  /+ *                      \|/+ *                       v1+ */+template <typename T, typename TNearPointLocator>+bool Triangulation<T, TNearPointLocator>::isFlipNeeded(+    const V2d<T>& v,+    const VertInd iV,+    const VertInd iV1,+    const VertInd iV2,+    const VertInd iV3) const+{+    const V2d<T>& v1 = vertices[iV1];+    const V2d<T>& v2 = vertices[iV2];+    const V2d<T>& v3 = vertices[iV3];+    if(m_superGeomType == SuperGeometryType::SuperTriangle)+    {+        // If flip-candidate edge touches super-triangle in-circumference+        // test has to be replaced with orient2d test against the line+        // formed by two non-artificial vertices (that don't belong to+        // super-triangle)+        if(iV < 3) // flip-candidate edge touches super-triangle+        {+            // does original edge also touch super-triangle?+            if(iV1 < 3)+                return locatePointLine(v1, v2, v3) ==+                       locatePointLine(v, v2, v3);+            if(iV3 < 3)+                return locatePointLine(v3, v1, v2) ==+                       locatePointLine(v, v1, v2);+            return false; // original edge does not touch super-triangle+        }+        if(iV2 < 3) // flip-candidate edge touches super-triangle+        {+            // does original edge also touch super-triangle?+            if(iV1 < 3)+                return locatePointLine(v1, v, v3) == locatePointLine(v2, v, v3);+            if(iV3 < 3)+                return locatePointLine(v3, v1, v) == locatePointLine(v2, v1, v);+            return false; // original edge does not touch super-triangle+        }+        // flip-candidate edge does not touch super-triangle+        if(iV1 < 3)+            return locatePointLine(v1, v2, v3) == locatePointLine(v, v2, v3);+        if(iV3 < 3)+            return locatePointLine(v3, v1, v2) == locatePointLine(v, v1, v2);+    }+    return isInCircumcircle(v, v1, v2, v3);+}++template <typename T, typename TNearPointLocator>+bool Triangulation<T, TNearPointLocator>::isFlipNeeded(+    const V2d<T>& v,+    const TriInd iT,+    const TriInd iTopo,+    const VertInd iV) const+{+    /*+     *                       v3         original edge: (v1, v3)+     *                      /|\   flip-candidate edge: (v,  v2)+     *                    /  |  \+     *                  /    |    \+     *                /      |      \+     * new vertex--> v       |       v2+     *                \      |      /+     *                  \    |    /+     *                    \  |  /+     *                      \|/+     *                       v1+     */+    const Triangle& tOpo = triangles[iTopo];+    const Index i = opposedVertexInd(tOpo, iT);+    const VertInd iV2 = tOpo.vertices[i];+    const VertInd iV1 = tOpo.vertices[cw(i)];+    const VertInd iV3 = tOpo.vertices[ccw(i)];++    // flip not needed if the original edge is fixed+    if(fixedEdges.count(Edge(iV1, iV3)))+        return false;++    return isFlipNeeded(v, iV, iV1, iV2, iV3);+}++/* Insert point into triangle: split into 3 triangles:+ *  - create 2 new triangles+ *  - re-use old triangle for the 3rd+ *                      v3+ *                    / | \+ *                   /  |  \ <-- original triangle (t)+ *                  /   |   \+ *              n3 /    |    \ n2+ *                /newT2|newT1\+ *               /      v      \+ *              /    __/ \__    \+ *             /  __/       \__  \+ *            / _/      t'     \_ \+ *          v1 ___________________ v2+ *                     n1+ */+template <typename T, typename TNearPointLocator>+std::stack<TriInd> Triangulation<T, TNearPointLocator>::insertPointInTriangle(+    const VertInd v,+    const TriInd iT)+{+    const TriInd iNewT1 = addTriangle();+    const TriInd iNewT2 = addTriangle();++    Triangle& t = triangles[iT];+    const array<VertInd, 3> vv = t.vertices;+    const array<TriInd, 3> nn = t.neighbors;+    const VertInd v1 = vv[0], v2 = vv[1], v3 = vv[2];+    const TriInd n1 = nn[0], n2 = nn[1], n3 = nn[2];+    // make two new triangles and convert current triangle to 3rd new+    // triangle+    using detail::arr3;+    triangles[iNewT1] = Triangle::make(arr3(v2, v3, v), arr3(n2, iNewT2, iT));+    triangles[iNewT2] = Triangle::make(arr3(v3, v1, v), arr3(n3, iT, iNewT1));+    t = Triangle::make(arr3(v1, v2, v), arr3(n1, iNewT1, iNewT2));+    // make and add a new vertex+    addAdjacentTriangles(v, iT, iNewT1, iNewT2);+    // adjust lists of adjacent triangles for v1, v2, v3+    addAdjacentTriangle(v1, iNewT2);+    addAdjacentTriangle(v2, iNewT1);+    removeAdjacentTriangle(v3, iT);+    addAdjacentTriangle(v3, iNewT1);+    addAdjacentTriangle(v3, iNewT2);+    // change triangle neighbor's neighbors to new triangles+    changeNeighbor(n2, iT, iNewT1);+    changeNeighbor(n3, iT, iNewT2);+    // return newly added triangles+    std::stack<TriInd> newTriangles;+    newTriangles.push(iT);+    newTriangles.push(iNewT1);+    newTriangles.push(iNewT2);+    return newTriangles;+}++/* Inserting a point on the edge between two triangles+ *    T1 (top)        v1+ *                   /|\+ *              n1 /  |  \ n4+ *               /    |    \+ *             /  T1' | Tnew1\+ *           v2-------v-------v4+ *             \ Tnew2| T2'  /+ *               \    |    /+ *              n2 \  |  / n3+ *                   \|/+ *   T2 (bottom)      v3+ */+template <typename T, typename TNearPointLocator>+std::stack<TriInd> Triangulation<T, TNearPointLocator>::insertPointOnEdge(+    const VertInd v,+    const TriInd iT1,+    const TriInd iT2)+{+    const TriInd iTnew1 = addTriangle();+    const TriInd iTnew2 = addTriangle();++    Triangle& t1 = triangles[iT1];+    Triangle& t2 = triangles[iT2];+    Index i = opposedVertexInd(t1, iT2);+    const VertInd v1 = t1.vertices[i];+    const VertInd v2 = t1.vertices[ccw(i)];+    const TriInd n1 = t1.neighbors[i];+    const TriInd n4 = t1.neighbors[cw(i)];+    i = opposedVertexInd(t2, iT1);+    const VertInd v3 = t2.vertices[i];+    const VertInd v4 = t2.vertices[ccw(i)];+    const TriInd n3 = t2.neighbors[i];+    const TriInd n2 = t2.neighbors[cw(i)];+    // add new triangles and change existing ones+    using detail::arr3;+    t1 = Triangle::make(arr3(v1, v2, v), arr3(n1, iTnew2, iTnew1));+    t2 = Triangle::make(arr3(v3, v4, v), arr3(n3, iTnew1, iTnew2));+    triangles[iTnew1] = Triangle::make(arr3(v1, v, v4), arr3(iT1, iT2, n4));+    triangles[iTnew2] = Triangle::make(arr3(v3, v, v2), arr3(iT2, iT1, n2));+    // make and add new vertex+    addAdjacentTriangles(v, iT1, iTnew2, iT2, iTnew1);+    // adjust neighboring triangles and vertices+    changeNeighbor(n4, iT1, iTnew1);+    changeNeighbor(n2, iT2, iTnew2);+    addAdjacentTriangle(v1, iTnew1);+    addAdjacentTriangle(v3, iTnew2);+    removeAdjacentTriangle(v2, iT2);+    addAdjacentTriangle(v2, iTnew2);+    removeAdjacentTriangle(v4, iT1);+    addAdjacentTriangle(v4, iTnew1);+    // return newly added triangles+    std::stack<TriInd> newTriangles;+    newTriangles.push(iT1);+    newTriangles.push(iTnew2);+    newTriangles.push(iT2);+    newTriangles.push(iTnew1);+    return newTriangles;+}++template <typename T, typename TNearPointLocator>+array<TriInd, 2>+Triangulation<T, TNearPointLocator>::trianglesAt(const V2d<T>& pos) const+{+    array<TriInd, 2> out = {noNeighbor, noNeighbor};+    for(TriInd i = TriInd(0); i < TriInd(triangles.size()); ++i)+    {+        const Triangle& t = triangles[i];+        const V2d<T>& v1 = vertices[t.vertices[0]];+        const V2d<T>& v2 = vertices[t.vertices[1]];+        const V2d<T>& v3 = vertices[t.vertices[2]];+        const PtTriLocation::Enum loc = locatePointTriangle(pos, v1, v2, v3);+        if(loc == PtTriLocation::Outside)+            continue;+        out[0] = i;+        if(isOnEdge(loc))+            out[1] = t.neighbors[edgeNeighbor(loc)];+        return out;+    }+    throw std::runtime_error("No triangle was found at position");+}++template <typename T, typename TNearPointLocator>+TriInd Triangulation<T, TNearPointLocator>::walkTriangles(+    const VertInd startVertex,+    const V2d<T>& pos) const+{+    // begin walk in search of triangle at pos+    TriInd currTri = vertTris[startVertex][0];+#ifdef CDT_USE_BOOST+    TriIndFlatUSet visited;+#else+    TriIndUSet visited;+#endif+    bool found = false;+    while(!found)+    {+        const Triangle& t = triangles[currTri];+        found = true;+        // stochastic offset to randomize which edge we check first+        const Index offset(detail::randGenerator() % 3);+        for(Index i_(0); i_ < Index(3); ++i_)+        {+            const Index i((i_ + offset) % 3);+            const V2d<T>& vStart = vertices[t.vertices[i]];+            const V2d<T>& vEnd = vertices[t.vertices[ccw(i)]];+            const PtLineLocation::Enum edgeCheck =+                locatePointLine(pos, vStart, vEnd);+            if(edgeCheck == PtLineLocation::Right &&+               t.neighbors[i] != noNeighbor &&+               visited.insert(t.neighbors[i]).second)+            {+                found = false;+                currTri = t.neighbors[i];+                break;+            }+        }+    }+    return currTri;+}++template <typename T, typename TNearPointLocator>+array<TriInd, 2> Triangulation<T, TNearPointLocator>::walkingSearchTrianglesAt(+    const V2d<T>& pos) const+{+    array<TriInd, 2> out = {noNeighbor, noNeighbor};+    // Query  for a vertex close to pos, to start the search+    const VertInd startVertex = m_nearPtLocator.nearPoint(pos, vertices);+    const TriInd iT = walkTriangles(startVertex, pos);+    // Finished walk, locate point in current triangle+    const Triangle& t = triangles[iT];+    const V2d<T>& v1 = vertices[t.vertices[0]];+    const V2d<T>& v2 = vertices[t.vertices[1]];+    const V2d<T>& v3 = vertices[t.vertices[2]];+    const PtTriLocation::Enum loc = locatePointTriangle(pos, v1, v2, v3);+    if(loc == PtTriLocation::Outside)+        throw std::runtime_error("No triangle was found at position");+    out[0] = iT;+    if(isOnEdge(loc))+        out[1] = t.neighbors[edgeNeighbor(loc)];+    return out;+}++/* Flip edge between T and Topo:+ *+ *                v4         | - old edge+ *               /|\         ~ - new edge+ *              / | \+ *          n3 /  T' \ n4+ *            /   |   \+ *           /    |    \+ *     T -> v1~~~~~~~~~v3 <- Topo+ *           \    |    /+ *            \   |   /+ *          n1 \Topo'/ n2+ *              \ | /+ *               \|/+ *                v2+ */+template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::flipEdge(+    const TriInd iT,+    const TriInd iTopo)+{+    Triangle& t = triangles[iT];+    Triangle& tOpo = triangles[iTopo];+    const array<TriInd, 3>& triNs = t.neighbors;+    const array<TriInd, 3>& triOpoNs = tOpo.neighbors;+    const array<VertInd, 3>& triVs = t.vertices;+    const array<VertInd, 3>& triOpoVs = tOpo.vertices;+    // find vertices and neighbors+    Index i = opposedVertexInd(t, iTopo);+    const VertInd v1 = triVs[i];+    const VertInd v2 = triVs[ccw(i)];+    const TriInd n1 = triNs[i];+    const TriInd n3 = triNs[cw(i)];+    i = opposedVertexInd(tOpo, iT);+    const VertInd v3 = triOpoVs[i];+    const VertInd v4 = triOpoVs[ccw(i)];+    const TriInd n4 = triOpoNs[i];+    const TriInd n2 = triOpoNs[cw(i)];+    // change vertices and neighbors+    using detail::arr3;+    t = Triangle::make(arr3(v4, v1, v3), arr3(n3, iTopo, n4));+    tOpo = Triangle::make(arr3(v2, v3, v1), arr3(n2, iT, n1));+    // adjust neighboring triangles and vertices+    changeNeighbor(n1, iT, iTopo);+    changeNeighbor(n4, iTopo, iT);+    // only adjust adjacent triangles if triangulation is not finalized:+    // can happen when called from outside on an already finalized triangulation+    if(!isFinalized())+    {+        addAdjacentTriangle(v1, iTopo);+        addAdjacentTriangle(v3, iT);+        removeAdjacentTriangle(v2, iT);+        removeAdjacentTriangle(v4, iTopo);+    }+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::changeNeighbor(+    const TriInd iT,+    const TriInd oldNeighbor,+    const TriInd newNeighbor)+{+    if(iT == noNeighbor)+        return;+    Triangle& t = triangles[iT];+    t.neighbors[neighborInd(t, oldNeighbor)] = newNeighbor;+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::addAdjacentTriangle(+    const VertInd iVertex,+    const TriInd iTriangle)+{+    vertTris[iVertex].push_back(iTriangle);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::addAdjacentTriangles(+    const VertInd iVertex,+    const TriInd iT1,+    const TriInd iT2,+    const TriInd iT3)+{+    TriIndVec& vTris = vertTris[iVertex];+    vTris.reserve(vTris.size() + 3);+    vTris.push_back(iT1);+    vTris.push_back(iT2);+    vTris.push_back(iT3);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::addAdjacentTriangles(+    const VertInd iVertex,+    const TriInd iT1,+    const TriInd iT2,+    const TriInd iT3,+    const TriInd iT4)+{+    TriIndVec& vTris = vertTris[iVertex];+    vTris.reserve(vTris.size() + 4);+    vTris.push_back(iT1);+    vTris.push_back(iT2);+    vTris.push_back(iT3);+    vTris.push_back(iT4);+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::removeAdjacentTriangle(+    const VertInd iVertex,+    const TriInd iTriangle)+{+    std::vector<TriInd>& tris = vertTris[iVertex];+    tris.erase(std::find(tris.begin(), tris.end(), iTriangle));+}++template <typename T, typename TNearPointLocator>+TriInd Triangulation<T, TNearPointLocator>::triangulatePseudopolygon(+    const VertInd ia,+    const VertInd ib,+    const std::vector<VertInd>::const_iterator pointsFirst,+    const std::vector<VertInd>::const_iterator pointsLast)+{+    if(pointsFirst == pointsLast)+        return pseudopolyOuterTriangle(ia, ib);+    // Find delaunay point+    const VertInd ic = findDelaunayPoint(ia, ib, pointsFirst, pointsLast);+    // Find pseudopolygons split by the delaunay point+    std::vector<VertInd>::const_iterator newLast = pointsFirst;+    while(*newLast != ic)+        ++newLast;+    const std::vector<VertInd>::const_iterator newFirst = newLast + 1;+    // triangulate splitted pseudo-polygons+    const TriInd iT2 = triangulatePseudopolygon(ic, ib, newFirst, pointsLast);+    const TriInd iT1 = triangulatePseudopolygon(ia, ic, pointsFirst, newLast);+    // add new triangle+    const Triangle t = {{ia, ib, ic}, {noNeighbor, iT2, iT1}};+    const TriInd iT = addTriangle(t);+    // adjust neighboring triangles and vertices+    if(iT1 != noNeighbor)+    {+        if(pointsFirst == newLast)+            changeNeighbor(iT1, ia, ic, iT);+        else+            triangles[iT1].neighbors[0] = iT;+    }+    if(iT2 != noNeighbor)+    {+        if(newFirst == pointsLast)+            changeNeighbor(iT2, ic, ib, iT);+        else+            triangles[iT2].neighbors[0] = iT;+    }+    addAdjacentTriangle(ia, iT);+    addAdjacentTriangle(ib, iT);+    addAdjacentTriangle(ic, iT);++    return iT;+}++template <typename T, typename TNearPointLocator>+VertInd Triangulation<T, TNearPointLocator>::findDelaunayPoint(+    const VertInd ia,+    const VertInd ib,+    const std::vector<VertInd>::const_iterator pointsFirst,+    const std::vector<VertInd>::const_iterator pointsLast) const+{+    assert(pointsFirst != pointsLast);+    const V2d<T>& a = vertices[ia];+    const V2d<T>& b = vertices[ib];+    VertInd ic = *pointsFirst;+    V2d<T> c = vertices[ic];+    typedef std::vector<VertInd>::const_iterator CIt;+    for(CIt it = pointsFirst + 1; it != pointsLast; ++it)+    {+        const V2d<T> v = vertices[*it];+        if(!isInCircumcircle(v, a, b, c))+            continue;+        ic = *it;+        c = vertices[ic];+    }+    return ic;+}++template <typename T, typename TNearPointLocator>+TriInd Triangulation<T, TNearPointLocator>::pseudopolyOuterTriangle(+    const VertInd ia,+    const VertInd ib) const+{+    const std::vector<TriInd>& aTris = vertTris[ia];+    const std::vector<TriInd>& bTris = vertTris[ib];+    typedef std::vector<TriInd>::const_iterator TriIndCit;+    for(TriIndCit it = aTris.begin(); it != aTris.end(); ++it)+        if(std::find(bTris.begin(), bTris.end(), *it) != bTris.end())+            return *it;+    return noNeighbor;+}++template <typename T, typename TNearPointLocator>+void Triangulation<T, TNearPointLocator>::insertVertices(+    const std::vector<V2d<T> >& newVertices)+{+    return insertVertices(+        newVertices.begin(), newVertices.end(), getX_V2d<T>, getY_V2d<T>);+}++template <typename T, typename TNearPointLocator>+bool Triangulation<T, TNearPointLocator>::isFinalized() const+{+    return vertTris.empty() && !vertices.empty();+}++template <typename T, typename TNearPointLocator>+unordered_map<TriInd, LayerDepth>+Triangulation<T, TNearPointLocator>::peelLayer(+    std::stack<TriInd> seeds,+    const LayerDepth layerDepth,+    std::vector<LayerDepth>& triDepths) const+{+    unordered_map<TriInd, LayerDepth> behindBoundary;+    while(!seeds.empty())+    {+        const TriInd iT = seeds.top();+        seeds.pop();+        triDepths[iT] = layerDepth;+        behindBoundary.erase(iT);+        const Triangle& t = triangles[iT];+        for(Index i(0); i < Index(3); ++i)+        {+            const Edge opEdge(t.vertices[ccw(i)], t.vertices[cw(i)]);+            const TriInd iN = t.neighbors[opoNbr(i)];+            if(iN == noNeighbor || triDepths[iN] <= layerDepth)+                continue;+            if(fixedEdges.count(opEdge))+            {+                const unordered_map<Edge, LayerDepth>::const_iterator cit =+                    overlapCount.find(opEdge);+                const LayerDepth triDepth = cit == overlapCount.end()+                                                ? layerDepth + 1+                                                : layerDepth + cit->second + 1;+                behindBoundary[iN] = triDepth;+                continue;+            }+            seeds.push(iN);+        }+    }+    return behindBoundary;+}++template <typename T, typename TNearPointLocator>+std::vector<LayerDepth>+Triangulation<T, TNearPointLocator>::calculateTriangleDepths() const+{+    std::vector<LayerDepth> triDepths(+        triangles.size(), std::numeric_limits<LayerDepth>::max());+    std::stack<TriInd> seeds(TriDeque(1, vertTris[0].front()));+    LayerDepth layerDepth = 0;+    LayerDepth deepestSeedDepth = 0;++    unordered_map<LayerDepth, TriIndUSet> seedsByDepth;+    do+    {+        const unordered_map<TriInd, LayerDepth>& newSeeds =+            peelLayer(seeds, layerDepth, triDepths);++        seedsByDepth.erase(layerDepth);+        typedef unordered_map<TriInd, LayerDepth>::const_iterator Iter;+        for(Iter it = newSeeds.begin(); it != newSeeds.end(); ++it)+        {+            deepestSeedDepth = std::max(deepestSeedDepth, it->second);+            seedsByDepth[it->second].insert(it->first);+        }+        const TriIndUSet& nextLayerSeeds = seedsByDepth[layerDepth + 1];+        seeds = std::stack<TriInd>(+            TriDeque(nextLayerSeeds.begin(), nextLayerSeeds.end()));+        ++layerDepth;+    } while(!seeds.empty() || deepestSeedDepth > layerDepth);++    return triDepths;+}++} // namespace CDT

@@ -0,0 +1,35 @@+typedef struct Vertex {+  double x;+  double y;+} VertexT;++typedef struct Edge {+  unsigned i;+  unsigned j;+} EdgeT;++typedef struct Triangle {+  unsigned i1;+  unsigned i2;+  unsigned i3;+} TriangleT;++typedef struct Triangulation {+  VertexT*   vertices;+  size_t     nvertices;+  TriangleT* triangles;+  size_t     ntriangles;+  EdgeT*     edges;+  size_t     nedges;+} TriangulationT;++typedef struct CTriangulation {+  VertexT*   vertices;+  size_t     nvertices;+  TriangleT* triangles;+  size_t     ntriangles;+  EdgeT*     edges;+  size_t     nedges;+  EdgeT*     fixededges;+  size_t     nfixededges;+} CTriangulationT;

@@ -0,0 +1,939 @@+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *+ *                                                                                 *+ * Copyright (c) 2019, William C. Lenthe                                           *+ * All rights reserved.                                                            *+ *                                                                                 *+ * Redistribution and use in source and binary forms, with or without              *+ * modification, are permitted provided that the following conditions are met:     *+ *                                                                                 *+ * 1. Redistributions of source code must retain the above copyright notice, this  *+ *    list of conditions and the following disclaimer.                             *+ *                                                                                 *+ * 2. 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.                       *+ *                                                                                 *+ * 3. Neither the name of the copyright holder 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.            *+ *                                                                                 *+ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */++#ifndef PREDICATES_H_INCLUDED+#define PREDICATES_H_INCLUDED++//@reference: https://www.cs.cmu.edu/~quake/robust.html++namespace  predicates {+	//@brief: geometric predicates using arbitrary precision arithmetic +	//@note : these are provided primarily for illustrative purposes and adaptive routines should be preferred+	namespace exact {+		//@brief   : determine if the 2d point c is above, on, or below the line defined by a and b+		//@param ax: X-coordinate of a+		//@param ay: Y-coordinate of a+		//@param bx: X-coordinate of b+		//@param by: Y-coordinate of b+		//@param cx: X-coordinate of c+		//@param cy: Y-coordinate of c+		//@return  : determinant of {{ax - cx, ay - cy}, {bx - cx, by - cy}}+		//@note    : positive, 0, negative result for c above, on, or below the line defined by a -> b+		template <typename T> T orient2d(T const ax, T const ay, T const bx, T const by, T const cx, T const cy);++		//@brief   : determine if the 2d point c is above, on, or below the line defined by a and b+		//@param pa: pointer to a as {x, y}+		//@param pb: pointer to b as {x, y}+		//@param pc: pointer to c as {x, y}+		//@return  : determinant of {{ax - cx, ay - cy}, {bx - cx, by - cy}}+		//@note    : positive, 0, negative result for c above, on, or below the line defined by a -> b+		template <typename T> T orient2d(T const*const pa, T const*const pb, T const*const pc);++		//@brief   : determine if the 2d point d is inside, on, or outside the circle defined by a, b, and c+		//@param ax: X-coordinate of a+		//@param ay: Y-coordinate of a+		//@param bx: X-coordinate of b+		//@param by: Y-coordinate of b+		//@param cx: X-coordinate of c+		//@param cy: Y-coordinate of c+		//@param dx: X-coordinate of d+		//@param dy: Y-coordinate of d+		//@return  : determinant of {{ax - dx, ay - dy, (ax - dx)^2 + (ay - dy)^2}, {bx - dx, by - dy, (bx - dx)^2 + (by - dy)^2}, {cx - dx, cy - dy, (cx - dx)^2 + (cy - dy)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T incircle(T const ax, T const ay, T const bx, T const by, T const cx, T const cy, T const dx, T const dy);++		//@brief   : determine if the 2d point d is inside, on, or outside the circle defined by a, b, and c+		//@param pa: pointer to a as {x, y}+		//@param pb: pointer to b as {x, y}+		//@param pc: pointer to c as {x, y}+		//@param pc: pointer to d as {x, y}+		//@return  : determinant of {{ax - dx, ay - dy, (ax - dx)^2 + (ay - dy)^2}, {bx - dx, by - dy, (bx - dx)^2 + (by - dy)^2}, {cx - dx, cy - dy, (cx - dx)^2 + (cy - dy)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T incircle(T const*const pa, T const*const pb, T const*const pc, T const*const pd);++		//@brief   : determine if the 3d point d is above, on, or below the plane defined by a, b, and c+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@return  : determinant of {{ax - dx, ay - dy, az - dz}, {bx - dx, by - dy, bz - dz}, {cx - dx, cy - dy, cz - dz}}+		//@note    : positive, 0, negative result for c above, on, or below the plane defined by a, b, and c+		template <typename T> T orient3d(T const*const pa, T const*const pb, T const*const pc, T const*const pd);+		+		//@brief   : determine if the 3d point e is inside, on, or outside the sphere defined by a, b, c, and d+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@param pe: pointer to e as {x, y, z}+		//@return  : determinant of {{ax - ex, ay - ey, az - ez, (ax - ex)^2 + (ay - ey)^2 + (az - ez)^2}, {bx - ex, by - ey, bz - ez, (bx - ex)^2 + (by - ey)^2 + (bz - ez)^2}, {cx - ex, cy - ey, cz - ez, (cx - ex)^2 + (cy - ey)^2 + (cz - ez)^2}, {dx - ex, dy - ey, dz - ez, (dx - ex)^2 + (dy - ey)^2 + (dz - ez)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T insphere(T const*const pa, T const*const pb, T const*const pc, T const*const pd, T const*const pe);+	}++	//@brief: geometric predicates using normal floating point arithmetic but falling back to arbitrary precision when needed+	//@note : these should have the same accuracy but are significantly faster when determinants are large+	namespace adaptive {+		//@brief   : determine if the 2d point c is above, on, or below the line defined by a and b+		//@param ax: X-coordinate of a+		//@param ay: Y-coordinate of a+		//@param bx: X-coordinate of b+		//@param by: Y-coordinate of b+		//@param cx: X-coordinate of c+		//@param cy: Y-coordinate of c+		//@return  : determinant of {{ax - cx, ay - cy}, {bx - cx, by - cy}}+		//@note    : positive, 0, negative result for c above, on, or below the line defined by a -> b+		template <typename T> T orient2d(T const ax, T const ay, T const bx, T const by, T const cx, T const cy);++		//@brief   : determine if the 2d point c is above, on, or below the line defined by a and b+		//@param pa: pointer to a as {x, y}+		//@param pb: pointer to b as {x, y}+		//@param pc: pointer to c as {x, y}+		//@return  : determinant of {{ax - cx, ay - cy}, {bx - cx, by - cy}}+		//@note    : positive, 0, negative result for c above, on, or below the line defined by a -> b+		template <typename T> T orient2d(T const*const pa, T const*const pb, T const*const pc);++		//@brief   : determine if the 2d point d is inside, on, or outside the circle defined by a, b, and c+		//@param ax: X-coordinate of a+		//@param ay: Y-coordinate of a+		//@param bx: X-coordinate of b+		//@param by: Y-coordinate of b+		//@param cx: X-coordinate of c+		//@param cy: Y-coordinate of c+		//@param dx: X-coordinate of d+		//@param dy: Y-coordinate of d+		//@return  : determinant of {{ax - dx, ay - dy, (ax - dx)^2 + (ay - dy)^2}, {bx - dx, by - dy, (bx - dx)^2 + (by - dy)^2}, {cx - dx, cy - dy, (cx - dx)^2 + (cy - dy)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T incircle(T const ax, T const ay, T const bx, T const by, T const cx, T const cy, T const dx, T const dy);++		//@brief   : determine if the 2d point d is inside, on, or outside the circle defined by a, b, and c+		//@param pa: pointer to a as {x, y}+		//@param pb: pointer to b as {x, y}+		//@param pc: pointer to c as {x, y}+		//@param pc: pointer to d as {x, y}+		//@return  : determinant of {{ax - dx, ay - dy, (ax - dx)^2 + (ay - dy)^2}, {bx - dx, by - dy, (bx - dx)^2 + (by - dy)^2}, {cx - dx, cy - dy, (cx - dx)^2 + (cy - dy)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T incircle(T const*const pa, T const*const pb, T const*const pc, T const*const pd);++		//@brief   : determine if the 3d point d is above, on, or below the plane defined by a, b, and c+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@return  : determinant of {{ax - dx, ay - dy, az - dz}, {bx - dx, by - dy, bz - dz}, {cx - dx, cy - dy, cz - dz}}+		//@note    : positive, 0, negative result for c above, on, or below the plane defined by a, b, and c+		template <typename T> T orient3d(T const*const pa, T const*const pb, T const*const pc, T const*const pd);++		//@brief   : determine if the 3d point e is inside, on, or outside the sphere defined by a, b, c, and d+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@param pe: pointer to e as {x, y, z}+		//@return  : determinant of {{ax - ex, ay - ey, az - ez, (ax - ex)^2 + (ay - ey)^2 + (az - ez)^2}, {bx - ex, by - ey, bz - ez, (bx - ex)^2 + (by - ey)^2 + (bz - ez)^2}, {cx - ex, cy - ey, cz - ez, (cx - ex)^2 + (cy - ey)^2 + (cz - ez)^2}, {dx - ex, dy - ey, dz - ez, (dx - ex)^2 + (dy - ey)^2 + (dz - ez)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T insphere(T const*const pa, T const*const pb, T const*const pc, T const*const pd, T const*const pe);+	}+}++#include <cmath>//abs, fma+#include <limits>+#include <utility>//pair+#include <numeric>//accumulate+#include <algorithm>//transform, copy_n, merge+#include <functional>//negate++// a macro based static assert for pre c++11+#define PREDICATES_PORTABLE_STATIC_ASSERT(condition, message) typedef char message[(condition) ? 1 : -1]++// check if c++11 is supported+#if !defined(__cplusplus) && !defined(_MSC_VER)+	PREDICATES_PORTABLE_STATIC_ASSERT(false, couldnt_parse_cxx_standard)+#endif+#if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)+	#define PREDICATES_CXX11_IS_SUPPORTED+#endif++// choose to use c++11 features or their backports+#ifdef PREDICATES_CXX11_IS_SUPPORTED+#include <array>+#include <type_traits>// is_same, enable_if+#undef PREDICATES_PORTABLE_STATIC_ASSERT+#define PREDICATES_TOKEN_TO_STRING1(x)  #x+#define PREDICATES_TOKEN_TO_STRING(x)  PREDICATES_TOKEN_TO_STRING1(x)+#define PREDICATES_PORTABLE_STATIC_ASSERT(condition, message) static_assert(condition, PREDICATES_TOKEN_TO_STRING(message))+namespace  predicates {+namespace stdx {+	using std::array;+	using std::copy_n;+}+#else+namespace  predicates {+namespace stdx {+	// array+	template<typename T, size_t N>+	class array {+		T buff[N];+	public:+		T& operator[](const size_t& i) { return buff[i]; }+		const T& operator[](const size_t& i) const { return buff[i]; }++		T       * data() { return buff; }+		T const * data() const { return buff; }++		T       *  begin() { return buff; }+		T const * cbegin() const { return buff; }+	};+	// copy_n+	template< class InputIt, class Size, class OutputIt>+	OutputIt copy_n(InputIt first, Size count, OutputIt result)+	{+		if (count > 0) {+			*result++ = *first;+			for (Size i = 1; i < count; ++i) {+				*result++ = *++first;+			}+		}+		return result;+	}+}+#endif // PREDICATES_CXX11_IS_SUPPORTED++namespace detail {+	template<typename T> class ExpansionBase;++	//@brief: class to exactly represent the result of a sequence of arithmetic operations as an sequence of values that sum to the result+	template<typename T, size_t N>+	class Expansion : private ExpansionBase<T>, public stdx::array<T, N> {+		private:+		public:+			size_t m_size;+			template <typename S> friend class ExpansionBase;//access for base class+			template <typename S, size_t M> friend class Expansion;//access for expansions of different size++			Expansion() : m_size(0) {}+			template <size_t M> Expansion& operator=(const Expansion<T, M>& e) {+				PREDICATES_PORTABLE_STATIC_ASSERT(M <= N, cannot_assign_a_larger_expansion_to_a_smaller_expansion);+				stdx::copy_n(e.cbegin(), e.size(), stdx::array<T, N>::begin());+				m_size = e.size();+				return *this;+			}++			//vector like convenience functions+			size_t size() const {return m_size;}+			bool empty() const {return 0 == m_size;}+			void push_back(const T v) {stdx::array<T, N>::operator[](m_size++) = v;}++		public:+			//estimates of expansion value+			T mostSignificant() const {return empty() ? T(0) : stdx::array<T, N>::operator[](m_size - 1);}+			T estimate() const {return std::accumulate(stdx::array<T, N>::cbegin(), stdx::array<T, N>::cbegin() + size(), T(0));}++			template <size_t M> Expansion<T, N+M> operator+(const Expansion<T, M>& f) const {+				Expansion<T, N+M> h;+				h.m_size = ExpansionBase<T>::ExpansionSum(this->data(), this->size(), f.data(), f.size(), h.data());+				return h;+			}++			void negate() {std::transform(stdx::array<T, N>::cbegin(), stdx::array<T, N>::cbegin() + size(), stdx::array<T, N>::begin(), std::negate<T>());}+			Expansion operator-() const {Expansion e = *this; e.negate(); return e;}+			template <size_t M> Expansion<T, N+M> operator-(const Expansion<T, M>& f) const {return operator+(-f);}++			Expansion<T, 2*N> operator*(const T b) const {+				Expansion<T, 2*N> h;+				h.m_size = ExpansionBase<T>::ScaleExpansion(this->data(), this->size(), b, h.data());+				return h;+			}+	};++	//std::fma is faster than dekker's product when the processor instruction is available+	#ifdef FP_FAST_FMAF+		static const bool fp_fast_fmaf = true;+	#else+		static const bool fp_fast_fmaf = false;+	#endif++	#ifdef FP_FAST_FMA+		static const bool fp_fast_fma = true;+	#else+		static const bool fp_fast_fma = false;+	#endif++	#ifdef FP_FAST_FMAL+		static const bool fp_fast_fmal = true;+	#else+		static const bool fp_fast_fmal = false;+	#endif++	#ifdef PREDICATES_CXX11_IS_SUPPORTED+	template <typename T> struct use_fma {static const bool value = (std::is_same<T, float>::value       && fp_fast_fmaf) ||+	                                                                (std::is_same<T, double>::value      && fp_fast_fma)  ||+	                                                                (std::is_same<T, long double>::value && fp_fast_fmal);};+	#endif++	//@brief  : helper function to sort by absolute value+	//@param a: lhs item to compare+	//@param b: rhs item to compare+	//@return : true if |a| < |b|+	//@note   : defined since lambda functions aren't allow in c++03+	template <typename T> bool absLess(const T& a, const T& b) {return std::abs(a) < std::abs(b);}++	template<typename T>+	class ExpansionBase {+		private:+			static const T Splitter;++			PREDICATES_PORTABLE_STATIC_ASSERT(std::numeric_limits<T>::is_iec559, Requires_IEC_559_IEEE_754_floating_point_type);+			PREDICATES_PORTABLE_STATIC_ASSERT(2 == std::numeric_limits<T>::radix, Requires_base_2_floating_point_type);++			//combine result + roundoff error into expansion+			static inline Expansion<T, 2> MakeExpansion(const T value, const T tail) {+				Expansion<T, 2> e;+				if(T(0) != tail) e.push_back(tail);+				if(T(0) != value) e.push_back(value);+				return e;+			}++		protected:+			//add 2 expansions+			static size_t ExpansionSum(T const * const e, const size_t n, T const * const f, const size_t m, T * const h) {+				std::merge(e, e + n, f, f + m, h, absLess<T>);+				if(m == 0) return n;+				if(n == 0) return m;+				size_t hIndex = 0;+				T Q = h[0];+				T Qnew = h[1] + Q;+				T hh = FastPlusTail(h[1], Q, Qnew);+				Q = Qnew;+				if(T(0) != hh) h[hIndex++] = hh;+				for(size_t g = 2; g != n + m; ++g) {+					Qnew = Q + h[g];+					hh = PlusTail(Q, h[g], Qnew);+					Q = Qnew;+					if(T(0) != hh) h[hIndex++] = hh;+				}+				if(T(0) != Q) h[hIndex++] = Q;+				return hIndex;+			}++			//scale an expansion by a constant+			static size_t ScaleExpansion(T const * const e, const size_t n, const T b, T * const h) {+				if(n == 0 || T(0) == b) return 0;+				size_t hIndex = 0;+				T Q = e[0] * b;+				const std::pair<T, T> bSplit = Split(b);+				T hh = MultTailPreSplit(e[0], b, bSplit, Q);+				if(T(0) != hh) h[hIndex++] = hh;+				for(size_t eIndex = 1; eIndex < n; ++eIndex) {+					T Ti = e[eIndex] * b;+					T ti = MultTailPreSplit(e[eIndex], b, bSplit, Ti);+					T Qi = Q + ti;+					hh = PlusTail(Q, ti, Qi);+					if(T(0) != hh) h[hIndex++] = hh;+					Q = Ti + Qi;+					hh = FastPlusTail(Ti, Qi, Q);+					if(T(0) != hh) h[hIndex++] = hh;+				}+				if(T(0) != Q) h[hIndex++] = Q;+				return hIndex;+			}+		+		public:+			//roundoff error of x = a + b+			static inline T PlusTail(const T a, const T b, const T x) {+				const T bVirtual = x - a;+				const T aVirtual = x - bVirtual;+				const T bRoundoff = b - bVirtual;+				const T aRoundoff = a - aVirtual;+				return aRoundoff + bRoundoff;+			}++			//roundoff error of x = a + b if |a| > |b|+			static inline T FastPlusTail(const T a, const T b, const T x) {+				const T bVirtual = x - a;+				return b - bVirtual;+			}++			//roundoff error of x = a - b+			static inline T MinusTail(const T a, const T b, const T x) {+				const T bVirtual = a - x;+				const T aVirtual = x + bVirtual;+				const T bRoundoff = bVirtual - b;+				const T aRoundoff = a - aVirtual;+				return aRoundoff + bRoundoff;+			}++			//split a into 2 nonoverlapping values+			static inline std::pair<T, T> Split(const T a) {+				const T c = a * Splitter;+				const T aBig = c - a;+				const T aHi = c - aBig;+				return std::pair<T, T>(aHi, a - aHi);+			}++			//roundoff error of x = a * b via dekkers product+			static inline T DekkersProduct(const T /*a*/, const std::pair<T, T> aSplit, const T /*b*/, const std::pair<T, T> bSplit, const T p) {+				T y = p - T(aSplit.first * bSplit.first);+				y -= T(aSplit.second * bSplit.first);+				y -= T(aSplit.first * bSplit.second);+				return T(aSplit.second * bSplit.second) - y;+			}++			//roundoff error of x = a * b+#ifdef PREDICATES_CXX11_IS_SUPPORTED+			template <typename S = T> static typename std::enable_if< use_fma<S>::value, S>::type MultTail(const T a, const T b, const T p) {return std::fma(a, b, -p);}+			template <typename S = T> static typename std::enable_if<!use_fma<S>::value, S>::type MultTail(const T a, const T b, const T p) {return DekkersProduct(a, Split(a), b, Split(b), p);}++			template <typename S = T> static typename std::enable_if< use_fma<S>::value, S>::type MultTailPreSplit(const T a, const T b, const std::pair<T, T> /*bSplit*/, const T p) {return std::fma(a, b, -p);}+			template <typename S = T> static typename std::enable_if<!use_fma<S>::value, S>::type MultTailPreSplit(const T a, const T b, const std::pair<T, T> bSplit, const T p) {return DekkersProduct(a, Split(a), b, bSplit, p);}+#else+			static T MultTail(const T a, const T b, const T p) {return DekkersProduct(a, Split(a), b, Split(b), p);}+			static T MultTailPreSplit(const T a, const T b, const std::pair<T, T> bSplit, const T p) {return DekkersProduct(a, Split(a), b, bSplit, p);}+#endif+			//expand a + b+			static inline Expansion<T, 2> Plus(const T a, const T b) {+				const T x = a + b;+				return MakeExpansion(x, PlusTail(a, b, x));+			}++			//expand a - b+			static inline Expansion<T, 2> Minus(const T a, const T b) {return Plus(a, -b);}++			//expand a * b+			static inline Expansion<T, 2> Mult(const T a, const T b) {+				const T x = a * b;+				return MakeExpansion(x, MultTail(a, b, x));+			}++			//expand the determinant of {{ax, ay}, {bx, by}} (unrolled Mult(ax, by) - Mult(ay, bx))+			static inline Expansion<T, 4> TwoTwoDiff(const T ax, const T by, const T ay, const T bx) {+				const T axby1 = ax * by;+				const T axby0 = MultTail(ax, by, axby1);+				const T bxay1 = bx * ay;+				const T bxay0 = MultTail(bx, ay, bxay1);+				const T _i0 = axby0 - bxay0;+				const T x0 = MinusTail(axby0, bxay0, _i0);+				const T _j = axby1 + _i0;+				const T _0 = PlusTail(axby1, _i0, _j);+				const T _i1 = _0 - bxay1;+				const T x1 = MinusTail(_0, bxay1, _i1);+				const T x3 = _j + _i1;+				const T x2 = PlusTail(_j, _i1, x3);+				Expansion<T, 4> e;+				if(T(0) != x0) e.push_back(x0);+				if(T(0) != x1) e.push_back(x1);+				if(T(0) != x2) e.push_back(x2);+				if(T(0) != x3) e.push_back(x3);+				return e;+			}++			//TwoTwoDiff checking for zeros to avoid extra splitting+			static inline Expansion<T, 4> TwoTwoDiffZeroCheck(const T ax, const T by, const T ay, const T bx) {+				Expansion<T, 4> e;+				if(T(0) == ax && T(0) == ay) return e;+				else if(T(0) == ax) e = Mult(ay, bx);+				else if(T(0) == ay) e = Mult(ax, by);+				else e = TwoTwoDiff(ax, by, ay, bx);+				return e;+			}++			//(a * b) * c checking for zeros+			static inline Expansion<T, 4> ThreeProd(const T a, const T b, const T c) {return (T(0) == a || T(0) == b || T(0) == c) ? Expansion<T, 4>() : Mult(a, b) * c;}+	};++	template <typename T> const T ExpansionBase<T>::Splitter = static_cast<T>(+#ifdef PREDICATES_CXX11_IS_SUPPORTED+		std::exp2((std::numeric_limits<T>::digits + std::numeric_limits<T>::digits%2)/2 + 1)+#else+		std::ldexp(T(1), (std::numeric_limits<T>::digits + std::numeric_limits<T>::digits%2)/2 + 1)+#endif+	);+}++	namespace exact {+		template <typename T> T orient2d(T const ax, T const ay, T const bx, T const by, T const cx, T const cy)+		{+			const detail::Expansion<T, 4> aterms = detail::ExpansionBase<T>::TwoTwoDiff(ax, by, ax, cy);+			const detail::Expansion<T, 4> bterms = detail::ExpansionBase<T>::TwoTwoDiff(bx, cy, bx, ay);+			const detail::Expansion<T, 4> cterms = detail::ExpansionBase<T>::TwoTwoDiff(cx, ay, cx, by);+			const detail::Expansion<T, 12> w = aterms + bterms + cterms;+			return w.mostSignificant();+		}++		template <typename T> T orient2d(T const*const pa, T const*const pb, T const*const pc) {+			return orient2d(pa[0], pa[1], pb[0], pb[1], pc[0], pc[1]);+		}++		template <typename T> T incircle(T const ax, T const ay, T const bx, T const by, T const cx, T const cy, T const dx, T const dy) {+			const detail::Expansion<T, 4> ab = detail::ExpansionBase<T>::TwoTwoDiff(ax, by, bx, ay);+			const detail::Expansion<T, 4> bc = detail::ExpansionBase<T>::TwoTwoDiff(bx, cy, cx, by);+			const detail::Expansion<T, 4> cd = detail::ExpansionBase<T>::TwoTwoDiff(cx, dy, dx, cy);+			const detail::Expansion<T, 4> da = detail::ExpansionBase<T>::TwoTwoDiff(dx, ay, ax, dy);+			const detail::Expansion<T, 4> ac = detail::ExpansionBase<T>::TwoTwoDiff(ax, cy, cx, ay);+			const detail::Expansion<T, 4> bd = detail::ExpansionBase<T>::TwoTwoDiff(bx, dy, dx, by);++			const detail::Expansion<T, 12> abc = ab + bc - ac;+			const detail::Expansion<T, 12> bcd = bc + cd - bd;+			const detail::Expansion<T, 12> cda = cd + da + ac;+			const detail::Expansion<T, 12> dab = da + ab + bd;++			const detail::Expansion<T, 96> adet = bcd * ax *  ax + bcd * ay *  ay;+			const detail::Expansion<T, 96> bdet = cda * bx * -bx + cda * by * -by;+			const detail::Expansion<T, 96> cdet = dab * cx *  cx + dab * cy *  cy;+			const detail::Expansion<T, 96> ddet = abc * dx * -dx + abc * dy * -dy;++			const detail::Expansion<T, 384> deter = (adet + bdet) + (cdet + ddet);+			return deter.mostSignificant();+		}++		template <typename T> T incircle(T const*const pa, T const*const pb, T const*const pc, T const*const pd) {+			return incircle(pa[0], pa[1], pb[0], pb[1], pc[0], pc[1], pd[0], pd[1]);+		}++		//@brief   : determine if the 3d point d is above, on, or below the plane defined by a, b, and c+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@return  : determinant of {{ax - dx, ay - dy, az - dz}, {bx - dx, by - dy, bz - dz}, {cx - dx, cy - dy, cz - dz}}+		//@note    : positive, 0, negative result for c above, on, or below the plane defined by a, b, and c+		template <typename T> T orient3d(T const*const pa, T const*const pb, T const*const pc, T const*const pd) {+			const detail::Expansion<T, 4> ab = detail::ExpansionBase<T>::TwoTwoDiff(pa[0], pb[1], pb[0], pa[1]);+			const detail::Expansion<T, 4> bc = detail::ExpansionBase<T>::TwoTwoDiff(pb[0], pc[1], pc[0], pb[1]);+			const detail::Expansion<T, 4> cd = detail::ExpansionBase<T>::TwoTwoDiff(pc[0], pd[1], pd[0], pc[1]);+			const detail::Expansion<T, 4> da = detail::ExpansionBase<T>::TwoTwoDiff(pd[0], pa[1], pa[0], pd[1]);+			const detail::Expansion<T, 4> ac = detail::ExpansionBase<T>::TwoTwoDiff(pa[0], pc[1], pc[0], pa[1]);+			const detail::Expansion<T, 4> bd = detail::ExpansionBase<T>::TwoTwoDiff(pb[0], pd[1], pd[0], pb[1]);++			const detail::Expansion<T, 12> abc = ab + bc - ac;+			const detail::Expansion<T, 12> bcd = bc + cd - bd;+			const detail::Expansion<T, 12> cda = cd + da + ac;+			const detail::Expansion<T, 12> dab = da + ab + bd;++			const detail::Expansion<T, 24> adet = bcd *  pa[2];+			const detail::Expansion<T, 24> bdet = cda * -pb[2];+			const detail::Expansion<T, 24> cdet = dab *  pc[2];+			const detail::Expansion<T, 24> ddet = abc * -pd[2];++			const detail::Expansion<T, 96> deter = (adet + bdet) + (cdet + ddet);+			return deter.mostSignificant();+		}++		//@brief   : determine if the 3d point e is inside, on, or outside the sphere defined by a, b, c, and d+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@param pe: pointer to e as {x, y, z}+		//@return  : determinant of {{ax - ex, ay - ey, az - ez, (ax - ex)^2 + (ay - ey)^2 + (az - ez)^2}, {bx - ex, by - ey, bz - ez, (bx - ex)^2 + (by - ey)^2 + (bz - ez)^2}, {cx - ex, cy - ey, cz - ez, (cx - ex)^2 + (cy - ey)^2 + (cz - ez)^2}, {dx - ex, dy - ey, dz - ez, (dx - ex)^2 + (dy - ey)^2 + (dz - ez)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T insphere(T const*const pa, T const*const pb, T const*const pc, T const*const pd, T const*const pe) {+			const detail::Expansion<T, 4> ab = detail::ExpansionBase<T>::TwoTwoDiff(pa[0], pb[1], pb[0], pa[1]);+			const detail::Expansion<T, 4> bc = detail::ExpansionBase<T>::TwoTwoDiff(pb[0], pc[1], pc[0], pb[1]);+			const detail::Expansion<T, 4> cd = detail::ExpansionBase<T>::TwoTwoDiff(pc[0], pd[1], pd[0], pc[1]);+			const detail::Expansion<T, 4> de = detail::ExpansionBase<T>::TwoTwoDiff(pd[0], pe[1], pe[0], pd[1]);+			const detail::Expansion<T, 4> ea = detail::ExpansionBase<T>::TwoTwoDiff(pe[0], pa[1], pa[0], pe[1]);+			const detail::Expansion<T, 4> ac = detail::ExpansionBase<T>::TwoTwoDiff(pa[0], pc[1], pc[0], pa[1]);+			const detail::Expansion<T, 4> bd = detail::ExpansionBase<T>::TwoTwoDiff(pb[0], pd[1], pd[0], pb[1]);+			const detail::Expansion<T, 4> ce = detail::ExpansionBase<T>::TwoTwoDiff(pc[0], pe[1], pe[0], pc[1]);+			const detail::Expansion<T, 4> da = detail::ExpansionBase<T>::TwoTwoDiff(pd[0], pa[1], pa[0], pd[1]);+			const detail::Expansion<T, 4> eb = detail::ExpansionBase<T>::TwoTwoDiff(pe[0], pb[1], pb[0], pe[1]);++			const detail::Expansion<T, 24> abc = bc * pa[2] + ac * -pb[2] + ab * pc[2];+			const detail::Expansion<T, 24> bcd = cd * pb[2] + bd * -pc[2] + bc * pd[2];+			const detail::Expansion<T, 24> cde = de * pc[2] + ce * -pd[2] + cd * pe[2];+			const detail::Expansion<T, 24> dea = ea * pd[2] + da * -pe[2] + de * pa[2];+			const detail::Expansion<T, 24> eab = ab * pe[2] + eb * -pa[2] + ea * pb[2];+			const detail::Expansion<T, 24> abd = bd * pa[2] + da *  pb[2] + ab * pd[2];+			const detail::Expansion<T, 24> bce = ce * pb[2] + eb *  pc[2] + bc * pe[2];+			const detail::Expansion<T, 24> cda = da * pc[2] + ac *  pd[2] + cd * pa[2];+			const detail::Expansion<T, 24> deb = eb * pd[2] + bd *  pe[2] + de * pb[2];+			const detail::Expansion<T, 24> eac = ac * pe[2] + ce *  pa[2] + ea * pc[2];++			const detail::Expansion<T, 96> bcde = (cde + bce) - (deb + bcd);+			const detail::Expansion<T, 96> cdea = (dea + cda) - (eac + cde);+			const detail::Expansion<T, 96> deab = (eab + deb) - (abd + dea);+			const detail::Expansion<T, 96> eabc = (abc + eac) - (bce + eab);+			const detail::Expansion<T, 96> abcd = (bcd + abd) - (cda + abc);++			const detail::Expansion<T, 1152> adet = bcde * pa[0] * pa[0] + bcde * pa[1] * pa[1] + bcde * pa[2] * pa[2];+			const detail::Expansion<T, 1152> bdet = cdea * pb[0] * pb[0] + cdea * pb[1] * pb[1] + cdea * pb[2] * pb[2];+			const detail::Expansion<T, 1152> cdet = deab * pc[0] * pc[0] + deab * pc[1] * pc[1] + deab * pc[2] * pc[2];+			const detail::Expansion<T, 1152> ddet = eabc * pd[0] * pd[0] + eabc * pd[1] * pd[1] + eabc * pd[2] * pd[2];+			const detail::Expansion<T, 1152> edet = abcd * pe[0] * pe[0] + abcd * pe[1] * pe[1] + abcd * pe[2] * pe[2];++			const detail::Expansion<T, 5760> deter = (adet + bdet) + ((cdet + ddet) + edet);+			return deter.mostSignificant();+		}+	}++	template <typename T>+	const T& Epsilon()+	{+		static const T epsilon = static_cast<T>(+#ifdef PREDICATES_CXX11_IS_SUPPORTED+			std::exp2(-std::numeric_limits<T>::digits)+#else+			std::ldexp(T(1), -std::numeric_limits<T>::digits)+#endif+		);+		return epsilon;+	}++	template <typename T>+	class Constants {+		public:+			static const T epsilon, resulterrbound;+			static const T ccwerrboundA, ccwerrboundB, ccwerrboundC;+			static const T o3derrboundA, o3derrboundB, o3derrboundC;+			static const T iccerrboundA, iccerrboundB, iccerrboundC;+			static const T isperrboundA, isperrboundB, isperrboundC;+	};++	template <typename T> const T Constants<T>::epsilon = Epsilon<T>();+	template <typename T> const T Constants<T>::resulterrbound = (T( 3) + T(   8) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::ccwerrboundA   = (T( 3) + T(  16) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::ccwerrboundB   = (T( 2) + T(  12) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::ccwerrboundC   = (T( 9) + T(  64) * Epsilon<T>()) * Epsilon<T>() * Epsilon<T>();+	template <typename T> const T Constants<T>::o3derrboundA   = (T( 7) + T(  56) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::o3derrboundB   = (T( 3) + T(  28) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::o3derrboundC   = (T(26) + T( 288) * Epsilon<T>()) * Epsilon<T>() * Epsilon<T>();+	template <typename T> const T Constants<T>::iccerrboundA   = (T(10) + T(  96) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::iccerrboundB   = (T( 4) + T(  48) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::iccerrboundC   = (T(44) + T( 576) * Epsilon<T>()) * Epsilon<T>() * Epsilon<T>();+	template <typename T> const T Constants<T>::isperrboundA   = (T(16) + T( 224) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::isperrboundB   = (T( 5) + T(  72) * Epsilon<T>()) * Epsilon<T>();+	template <typename T> const T Constants<T>::isperrboundC   = (T(71) + T(1408) * Epsilon<T>()) * Epsilon<T>() * Epsilon<T>();++	namespace adaptive {+		template <typename T> T orient2d(T const ax, T const ay, T const bx, T const by, T const cx, T const cy) {+			const T acx = ax - cx;+			const T bcx = bx - cx;+			const T acy = ay - cy;+			const T bcy = by - cy;+			const T detleft = acx * bcy;+			const T detright = acy * bcx;+			T det = detleft - detright;+			if((detleft < 0) != (detright < 0)) return det;+			if(T(0) == detleft || T(0) == detright) return det;++			const T detsum = std::abs(detleft + detright);+			T errbound = Constants<T>::ccwerrboundA * detsum;+			if(std::abs(det) >= std::abs(errbound)) return det;++			const detail::Expansion<T, 4> B = detail::ExpansionBase<T>::TwoTwoDiff(acx, bcy, acy, bcx);+			det = B.estimate();+			errbound = Constants<T>::ccwerrboundB * detsum;+			if(std::abs(det) >= std::abs(errbound)) return det;++			const T acxtail = detail::ExpansionBase<T>::MinusTail(ax, cx, acx);+			const T bcxtail = detail::ExpansionBase<T>::MinusTail(bx, cx, bcx);+			const T acytail = detail::ExpansionBase<T>::MinusTail(ay, cy, acy);+			const T bcytail = detail::ExpansionBase<T>::MinusTail(by, cy, bcy);+			if(T(0) == acxtail && T(0) == bcxtail && T(0) == acytail && T(0) == bcytail) return det;++			errbound = Constants<T>::ccwerrboundC * detsum + Constants<T>::resulterrbound * std::abs(det);+			det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail);+			if(std::abs(det) >= std::abs(errbound)) return det;++			const detail::Expansion<T, 16> D = ((B + detail::ExpansionBase<T>::TwoTwoDiff(acxtail, bcy, acytail, bcx)) + detail::ExpansionBase<T>::TwoTwoDiff(acx, bcytail, acy, bcxtail)) + detail::ExpansionBase<T>::TwoTwoDiff(acxtail, bcytail, acytail, bcxtail);+			return D.mostSignificant();+		}++		template <typename T> T orient2d(T const*const pa, T const*const pb, T const*const pc) {+			return orient2d(pa[0], pa[1], pb[0], pb[1], pc[0], pc[1]);+		}++		template <typename T> T incircle(T const ax, T const ay, T const bx, T const by, T const cx, T const cy, T const dx, T const dy) {+			const T adx = ax - dx;+			const T bdx = bx - dx;+			const T cdx = cx - dx;+			const T ady = ay - dy;+			const T bdy = by - dy;+			const T cdy = cy - dy;+			const T bdxcdy = bdx * cdy;+			const T cdxbdy = cdx * bdy;+			const T cdxady = cdx * ady;+			const T adxcdy = adx * cdy;+			const T adxbdy = adx * bdy;+			const T bdxady = bdx * ady;+			const T alift = adx * adx + ady * ady;+			const T blift = bdx * bdx + bdy * bdy;+			const T clift = cdx * cdx + cdy * cdy;+			T det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);+			const T permanent = (std::abs(bdxcdy) + std::abs(cdxbdy)) * alift+			                  + (std::abs(cdxady) + std::abs(adxcdy)) * blift+			                  + (std::abs(adxbdy) + std::abs(bdxady)) * clift;+			T errbound = Constants<T>::iccerrboundA * permanent;+			if(std::abs(det) >= std::abs(errbound)) return det;++			const detail::Expansion<T, 4> bc = detail::ExpansionBase<T>::TwoTwoDiff(bdx, cdy, cdx, bdy);+			const detail::Expansion<T, 4> ca = detail::ExpansionBase<T>::TwoTwoDiff(cdx, ady, adx, cdy);+			const detail::Expansion<T, 4> ab = detail::ExpansionBase<T>::TwoTwoDiff(adx, bdy, bdx, ady);+			const detail::Expansion<T, 32> adet = bc * adx * adx + bc * ady * ady;+			const detail::Expansion<T, 32> bdet = ca * bdx * bdx + ca * bdy * bdy;+			const detail::Expansion<T, 32> cdet = ab * cdx * cdx + ab * cdy * cdy;+			const detail::Expansion<T, 96> fin1 = adet + bdet + cdet;+			det = fin1.estimate();+			errbound = Constants<T>::iccerrboundB * permanent;+			if(std::abs(det) >= std::abs(errbound)) return det;++			const T adxtail = detail::ExpansionBase<T>::MinusTail(ax, dx, adx);+			const T adytail = detail::ExpansionBase<T>::MinusTail(ay, dy, ady);+			const T bdxtail = detail::ExpansionBase<T>::MinusTail(bx, dx, bdx);+			const T bdytail = detail::ExpansionBase<T>::MinusTail(by, dy, bdy);+			const T cdxtail = detail::ExpansionBase<T>::MinusTail(cx, dx, cdx);+			const T cdytail = detail::ExpansionBase<T>::MinusTail(cy, dy, cdy);+			if(T(0) == adxtail && T(0) == bdxtail && T(0) == cdxtail && T(0) == adytail && T(0) == bdytail && T(0) == cdytail) return det;++			errbound = Constants<T>::iccerrboundC * permanent + Constants<T>::resulterrbound * std::abs(det);+			det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail))+			    +   (bdx * cdy - bdy * cdx) *  (adx * adxtail + ady * adytail) * T(2))+			    +  ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail))+			    +   (cdx * ady - cdy * adx) *  (bdx * bdxtail + bdy * bdytail) * T(2))+			    +  ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail))+			    +   (adx * bdy - ady * bdx) *  (cdx * cdxtail + cdy * cdytail) * T(2));+			if(std::abs(det) >= std::abs(errbound)) return det;+			return exact::incircle(ax, ay, bx, by, cx, cy, dx, dy);+		}++		template <typename T> T incircle(T const*const pa, T const*const pb, T const*const pc, T const*const pd) {+			return incircle(pa[0], pa[1], pb[0], pb[1], pc[0], pc[1], pd[0], pd[1]);+		}++		//@brief   : determine if the 3d point d is above, on, or below the plane defined by a, b, and c+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@return  : determinant of {{ax - dx, ay - dy, az - dz}, {bx - dx, by - dy, bz - dz}, {cx - dx, cy - dy, cz - dz}}+		//@note    : positive, 0, negative result for c above, on, or below the plane defined by a, b, and c+		template <typename T> T orient3d(T const*const pa, T const*const pb, T const*const pc, T const*const pd) {+			const T adx = pa[0] - pd[0];+			const T bdx = pb[0] - pd[0];+			const T cdx = pc[0] - pd[0];+			const T ady = pa[1] - pd[1];+			const T bdy = pb[1] - pd[1];+			const T cdy = pc[1] - pd[1];+			const T adz = pa[2] - pd[2];+			const T bdz = pb[2] - pd[2];+			const T cdz = pc[2] - pd[2];+			const T bdxcdy = bdx * cdy;+			const T cdxbdy = cdx * bdy;+			const T cdxady = cdx * ady;+			const T adxcdy = adx * cdy;+			const T adxbdy = adx * bdy;+			const T bdxady = bdx * ady;+			T det = adz * (bdxcdy - cdxbdy) + bdz * (cdxady - adxcdy) + cdz * (adxbdy - bdxady);+			const T permanent = (std::abs(bdxcdy) + std::abs(cdxbdy)) * std::abs(adz) + (std::abs(cdxady) + std::abs(adxcdy)) * std::abs(bdz) + (std::abs(adxbdy) + std::abs(bdxady)) * std::abs(cdz);+			T errbound = Constants<T>::o3derrboundA * permanent;+			if(std::abs(det) >= std::abs(errbound)) return det;++			const detail::Expansion<T, 4> bc = detail::ExpansionBase<T>::TwoTwoDiff(bdx, cdy, cdx, bdy);+			const detail::Expansion<T, 4> ca = detail::ExpansionBase<T>::TwoTwoDiff(cdx, ady, adx, cdy);+			const detail::Expansion<T, 4> ab = detail::ExpansionBase<T>::TwoTwoDiff(adx, bdy, bdx, ady);+			const detail::Expansion<T, 24> fin1 = (bc * adz + ca * bdz) + ab * cdz;+			det = fin1.estimate();+			errbound = Constants<T>::o3derrboundB * permanent;+			if(std::abs(det) >= std::abs(errbound)) return det;++			const T adxtail = detail::ExpansionBase<T>::MinusTail(pa[0], pd[0], adx);+			const T bdxtail = detail::ExpansionBase<T>::MinusTail(pb[0], pd[0], bdx);+			const T cdxtail = detail::ExpansionBase<T>::MinusTail(pc[0], pd[0], cdx);+			const T adytail = detail::ExpansionBase<T>::MinusTail(pa[1], pd[1], ady);+			const T bdytail = detail::ExpansionBase<T>::MinusTail(pb[1], pd[1], bdy);+			const T cdytail = detail::ExpansionBase<T>::MinusTail(pc[1], pd[1], cdy);+			const T adztail = detail::ExpansionBase<T>::MinusTail(pa[2], pd[2], adz);+			const T bdztail = detail::ExpansionBase<T>::MinusTail(pb[2], pd[2], bdz);+			const T cdztail = detail::ExpansionBase<T>::MinusTail(pc[2], pd[2], cdz);+			if(T(0) == adxtail && T(0) == adytail && T(0) == adztail &&+			   T(0) == bdxtail && T(0) == bdytail && T(0) == bdztail &&+			   T(0) == cdxtail && T(0) == cdytail && T(0) == cdztail) return det;++			errbound = Constants<T>::o3derrboundC * permanent + Constants<T>::resulterrbound * std::abs(det);+			det += (adz * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail)) + adztail * (bdx * cdy - bdy * cdx))+			    +  (bdz * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail)) + bdztail * (cdx * ady - cdy * adx))+			    +  (cdz * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail)) + cdztail * (adx * bdy - ady * bdx));+			if(std::abs(det) >= std::abs(errbound)) return det;++			const detail::Expansion<T, 8> bct = detail::ExpansionBase<T>::TwoTwoDiffZeroCheck(bdxtail, cdy, bdytail, cdx) + detail::ExpansionBase<T>::TwoTwoDiffZeroCheck(cdytail, bdx, cdxtail, bdy);+			const detail::Expansion<T, 8> cat = detail::ExpansionBase<T>::TwoTwoDiffZeroCheck(cdxtail, ady, cdytail, adx) + detail::ExpansionBase<T>::TwoTwoDiffZeroCheck(adytail, cdx, adxtail, cdy);+			const detail::Expansion<T, 8> abt = detail::ExpansionBase<T>::TwoTwoDiffZeroCheck(adxtail, bdy, adytail, bdx) + detail::ExpansionBase<T>::TwoTwoDiffZeroCheck(bdytail, adx, bdxtail, ady);+			const detail::Expansion<T, 192> fin2 = fin1 + bct * adz + cat * bdz + abt * cdz + bc * adztail + ca * bdztail + ab * cdztail+			                                     + detail::ExpansionBase<T>::ThreeProd( adxtail, bdytail, cdz) + detail::ExpansionBase<T>::ThreeProd( adxtail, bdytail, cdztail)+			                                     + detail::ExpansionBase<T>::ThreeProd(-adxtail, cdytail, bdz) + detail::ExpansionBase<T>::ThreeProd(-adxtail, cdytail, bdztail)+			                                     + detail::ExpansionBase<T>::ThreeProd( bdxtail, cdytail, adz) + detail::ExpansionBase<T>::ThreeProd( bdxtail, cdytail, adztail)+			                                     + detail::ExpansionBase<T>::ThreeProd(-bdxtail, adytail, cdz) + detail::ExpansionBase<T>::ThreeProd(-bdxtail, adytail, cdztail)+			                                     + detail::ExpansionBase<T>::ThreeProd( cdxtail, adytail, bdz) + detail::ExpansionBase<T>::ThreeProd( cdxtail, adytail, bdztail)+			                                     + detail::ExpansionBase<T>::ThreeProd(-cdxtail, bdytail, adz) + detail::ExpansionBase<T>::ThreeProd(-cdxtail, bdytail, adztail)+			                                     + bct * adztail + cat * bdztail + abt * cdztail;+			return fin2.mostSignificant();+		}++		//@brief   : determine if the 3d point e is inside, on, or outside the sphere defined by a, b, c, and d+		//@param pa: pointer to a as {x, y, z}+		//@param pb: pointer to b as {x, y, z}+		//@param pc: pointer to c as {x, y, z}+		//@param pd: pointer to d as {x, y, z}+		//@param pe: pointer to e as {x, y, z}+		//@return  : determinant of {{ax - ex, ay - ey, az - ez, (ax - ex)^2 + (ay - ey)^2 + (az - ez)^2}, {bx - ex, by - ey, bz - ez, (bx - ex)^2 + (by - ey)^2 + (bz - ez)^2}, {cx - ex, cy - ey, cz - ez, (cx - ex)^2 + (cy - ey)^2 + (cz - ez)^2}, {dx - ex, dy - ey, dz - ez, (dx - ex)^2 + (dy - ey)^2 + (dz - ez)^2}}+		//@note    : positive, 0, negative result for d inside, on, or outside the circle defined by a, b, and c+		template <typename T> T insphere(T const*const pa, T const*const pb, T const*const pc, T const*const pd, T const*const pe) {+			T permanent;+			const T aex = pa[0] - pe[0];+			const T bex = pb[0] - pe[0];+			const T cex = pc[0] - pe[0];+			const T dex = pd[0] - pe[0];+			const T aey = pa[1] - pe[1];+			const T bey = pb[1] - pe[1];+			const T cey = pc[1] - pe[1];+			const T dey = pd[1] - pe[1];+			const T aez = pa[2] - pe[2];+			const T bez = pb[2] - pe[2];+			const T cez = pc[2] - pe[2];+			const T dez = pd[2] - pe[2];+			{+				const T aexbey = aex * bey;+				const T bexaey = bex * aey;+				const T bexcey = bex * cey;+				const T cexbey = cex * bey;+				const T cexdey = cex * dey;+				const T dexcey = dex * cey;+				const T dexaey = dex * aey;+				const T aexdey = aex * dey;+				const T aexcey = aex * cey;+				const T cexaey = cex * aey;+				const T bexdey = bex * dey;+				const T dexbey = dex * bey;+				const T ab = aexbey - bexaey;+				const T bc = bexcey - cexbey;+				const T cd = cexdey - dexcey;+				const T da = dexaey - aexdey;+				const T ac = aexcey - cexaey;+				const T bd = bexdey - dexbey;+				const T abc = aez * bc - bez * ac + cez * ab;+				const T bcd = bez * cd - cez * bd + dez * bc;+				const T cda = cez * da + dez * ac + aez * cd;+				const T dab = dez * ab + aez * bd + bez * da;+				const T alift = aex * aex + aey * aey + aez * aez;+				const T blift = bex * bex + bey * bey + bez * bez;+				const T clift = cex * cex + cey * cey + cez * cez;+				const T dlift = dex * dex + dey * dey + dez * dez;+				const T det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);+				const T aezplus = std::abs(aez);+				const T bezplus = std::abs(bez);+				const T cezplus = std::abs(cez);+				const T dezplus = std::abs(dez);+				const T aexbeyplus = std::abs(aexbey);+				const T bexaeyplus = std::abs(bexaey);+				const T bexceyplus = std::abs(bexcey);+				const T cexbeyplus = std::abs(cexbey);+				const T cexdeyplus = std::abs(cexdey);+				const T dexceyplus = std::abs(dexcey);+				const T dexaeyplus = std::abs(dexaey);+				const T aexdeyplus = std::abs(aexdey);+				const T aexceyplus = std::abs(aexcey);+				const T cexaeyplus = std::abs(cexaey);+				const T bexdeyplus = std::abs(bexdey);+				const T dexbeyplus = std::abs(dexbey);+				permanent = ((cexdeyplus + dexceyplus) * bezplus + (dexbeyplus + bexdeyplus) * cezplus + (bexceyplus + cexbeyplus) * dezplus) * alift+				          + ((dexaeyplus + aexdeyplus) * cezplus + (aexceyplus + cexaeyplus) * dezplus + (cexdeyplus + dexceyplus) * aezplus) * blift+				          + ((aexbeyplus + bexaeyplus) * dezplus + (bexdeyplus + dexbeyplus) * aezplus + (dexaeyplus + aexdeyplus) * bezplus) * clift+				          + ((bexceyplus + cexbeyplus) * aezplus + (cexaeyplus + aexceyplus) * bezplus + (aexbeyplus + bexaeyplus) * cezplus) * dlift;+				const T errbound = Constants<T>::isperrboundA * permanent;+				if(std::abs(det) >= std::abs(errbound)) return det;+			}++			const detail::Expansion<T, 4> ab = detail::ExpansionBase<T>::TwoTwoDiff(aex, bey, bex, aey);+			const detail::Expansion<T, 4> bc = detail::ExpansionBase<T>::TwoTwoDiff(bex, cey, cex, bey);+			const detail::Expansion<T, 4> cd = detail::ExpansionBase<T>::TwoTwoDiff(cex, dey, dex, cey);+			const detail::Expansion<T, 4> da = detail::ExpansionBase<T>::TwoTwoDiff(dex, aey, aex, dey);+			const detail::Expansion<T, 4> ac = detail::ExpansionBase<T>::TwoTwoDiff(aex, cey, cex, aey);+			const detail::Expansion<T, 4> bd = detail::ExpansionBase<T>::TwoTwoDiff(bex, dey, dex, bey);+			const detail::Expansion<T, 24> temp24a = bc * dez + (cd * bez + bd * -cez);+			const detail::Expansion<T, 24> temp24b = cd * aez + (da * cez + ac *  dez);+			const detail::Expansion<T, 24> temp24c = da * bez + (ab * dez + bd *  aez);+			const detail::Expansion<T, 24> temp24d = ab * cez + (bc * aez + ac * -bez);+			const detail::Expansion<T, 288> adet = temp24a * aex * -aex + temp24a * aey * -aey + temp24a * aez * -aez;+			const detail::Expansion<T, 288> bdet = temp24b * bex *  bex + temp24b * bey *  bey + temp24b * bez *  bez;+			const detail::Expansion<T, 288> cdet = temp24c * cex * -cex + temp24c * cey * -cey + temp24c * cez * -cez;+			const detail::Expansion<T, 288> ddet = temp24d * dex *  dex + temp24d * dey *  dey + temp24d * dez *  dez;+			const detail::Expansion<T, 1152> fin1 = (adet + bdet) + (cdet + ddet);+			T det = fin1.estimate();+			T errbound = Constants<T>::isperrboundB * permanent;+			if(std::abs(det) >= std::abs(errbound)) return det;++			const T aextail = detail::ExpansionBase<T>::MinusTail(pa[0], pe[0], aex);+			const T aeytail = detail::ExpansionBase<T>::MinusTail(pa[1], pe[1], aey);+			const T aeztail = detail::ExpansionBase<T>::MinusTail(pa[2], pe[2], aez);+			const T bextail = detail::ExpansionBase<T>::MinusTail(pb[0], pe[0], bex);+			const T beytail = detail::ExpansionBase<T>::MinusTail(pb[1], pe[1], bey);+			const T beztail = detail::ExpansionBase<T>::MinusTail(pb[2], pe[2], bez);+			const T cextail = detail::ExpansionBase<T>::MinusTail(pc[0], pe[0], cex);+			const T ceytail = detail::ExpansionBase<T>::MinusTail(pc[1], pe[1], cey);+			const T ceztail = detail::ExpansionBase<T>::MinusTail(pc[2], pe[2], cez);+			const T dextail = detail::ExpansionBase<T>::MinusTail(pd[0], pe[0], dex);+			const T deytail = detail::ExpansionBase<T>::MinusTail(pd[1], pe[1], dey);+			const T deztail = detail::ExpansionBase<T>::MinusTail(pd[2], pe[2], dez);+			if (T(0) == aextail && T(0) == aeytail && T(0) == aeztail &&+			    T(0) == bextail && T(0) == beytail && T(0) == beztail &&+			    T(0) == cextail && T(0) == ceytail && T(0) == ceztail &&+			    T(0) == dextail && T(0) == deytail && T(0) == deztail) return det;++			errbound = Constants<T>::isperrboundC * permanent + Constants<T>::resulterrbound * std::abs(det);+			const T abeps = (aex * beytail + bey * aextail) - (aey * bextail + bex * aeytail);+			const T bceps = (bex * ceytail + cey * bextail) - (bey * cextail + cex * beytail);+			const T cdeps = (cex * deytail + dey * cextail) - (cey * dextail + dex * ceytail);+			const T daeps = (dex * aeytail + aey * dextail) - (dey * aextail + aex * deytail);+			const T aceps = (aex * ceytail + cey * aextail) - (aey * cextail + cex * aeytail);+			const T bdeps = (bex * deytail + dey * bextail) - (bey * dextail + dex * beytail);+			const T ab3 = ab.mostSignificant();+			const T bc3 = bc.mostSignificant();+			const T cd3 = cd.mostSignificant();+			const T da3 = da.mostSignificant();+			const T ac3 = ac.mostSignificant();+			const T bd3 = bd.mostSignificant();+			det += ( ( (bex * bex + bey * bey + bez * bez) * ((cez * daeps + dez * aceps + aez * cdeps) + (ceztail * da3 + deztail * ac3 + aeztail * cd3))+			         + (dex * dex + dey * dey + dez * dez) * ((aez * bceps - bez * aceps + cez * abeps) + (aeztail * bc3 - beztail * ac3 + ceztail * ab3)) )+			       - ( (aex * aex + aey * aey + aez * aez) * ((bez * cdeps - cez * bdeps + dez * bceps) + (beztail * cd3 - ceztail * bd3 + deztail * bc3))+			         + (cex * cex + cey * cey + cez * cez) * ((dez * abeps + aez * bdeps + bez * daeps) + (deztail * ab3 + aeztail * bd3 + beztail * da3)) ) )+			    + T(2) * ( ( (bex * bextail + bey * beytail + bez * beztail) * (cez * da3 + dez * ac3 + aez * cd3)+			               + (dex * dextail + dey * deytail + dez * deztail) * (aez * bc3 - bez * ac3 + cez * ab3))+			             - ( (aex * aextail + aey * aeytail + aez * aeztail) * (bez * cd3 - cez * bd3 + dez * bc3)+			               + (cex * cextail + cey * ceytail + cez * ceztail) * (dez * ab3 + aez * bd3 + bez * da3)));+			if(std::abs(det) >= std::abs(errbound)) return det;+			return exact::insphere(pa, pb, pc, pd, pe);+		}+	}+}++#endif

@@ -0,0 +1,55 @@+#ifndef REMOVE_AT_HPP+#define REMOVE_AT_HPP++// check if c++11 is supported+#if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)+#define REMOVE_AT_CXX11_IS_SUPPORTED+#elif !defined(__cplusplus) && !defined(_MSC_VER)+typedef char couldnt_parse_cxx_standard[-1];+#endif++#include <algorithm>+#include <iterator>++/*!+ * Remove elements in the range [first; last) with indices from the sorted+ * unique range [ii_first, ii_last)+ */+template <class ForwardIt, class SortUniqIndsFwdIt>+inline ForwardIt remove_at(+    ForwardIt first,+    ForwardIt last,+    SortUniqIndsFwdIt ii_first,+    SortUniqIndsFwdIt ii_last)+{+    if(ii_first == ii_last) // no indices-to-remove are given+        return last;+    typedef typename std::iterator_traits<ForwardIt>::difference_type diff_t;+    typedef typename std::iterator_traits<SortUniqIndsFwdIt>::value_type ind_t;+    ForwardIt destination = first + static_cast<diff_t>(*ii_first);+    while(ii_first != ii_last)+    {+        // advance to an index after a chunk of elements-to-keep+        for(ind_t cur = *ii_first++; ii_first != ii_last; ++ii_first)+        {+            const ind_t nxt = *ii_first;+            if(nxt - cur > 1)+                break;+            cur = nxt;+        }+        // move the chunk of elements-to-keep to new destination+        const ForwardIt source_first =+            first + static_cast<diff_t>(*(ii_first - 1)) + 1;+        const ForwardIt source_last =+            ii_first != ii_last ? first + static_cast<diff_t>(*ii_first) : last;+#ifdef REMOVE_AT_CXX11_IS_SUPPORTED+        std::move(source_first, source_last, destination);+#else+        std::copy(source_first, source_last, destination); // c++98 version+#endif+        destination += source_last - source_first;+    }+    return destination;+}++#endif // REMOVE_AT_HPP

@@ -1,6 +1,6 @@ cabal-version:       2.2 name:                hcdt-version:             0.1.0.2+version:             0.1.0.3 synopsis:            2d Delaunay triangulation description:              This library performs the constrained or unconstrained 2d Delaunay triangulation.@@ -31,17 +31,17 @@                      , containers >= 0.6.4.1 && < 0.7                      , indexed-traversable >= 0.1.2 && < 0.2   include-dirs:        cpp-  includes:            ./cpp/CDT.h-                     , ./cpp/CDT.hpp-                     , ./cpp/CDTUtils.h-                     , ./cpp/CDTUtils.hpp-                     , ./cpp/hcdt.hpp-                     , ./cpp/KDTree.h-                     , ./cpp/LocatorKDTree.h-                     , ./cpp/predicates.h-                     , ./cpp/remove_at.hpp-                     , ./cpp/Triangulation.h-                     , ./cpp/Triangulation.hpp+  install-includes:    cpp/CDT.h+                     , cpp/CDT.hpp+                     , cpp/CDTUtils.h+                     , cpp/CDTUtils.hpp+                     , cpp/hcdt.hpp+                     , cpp/KDTree.h+                     , cpp/LocatorKDTree.h+                     , cpp/predicates.h+                     , cpp/remove_at.hpp+                     , cpp/Triangulation.h+                     , cpp/Triangulation.hpp   C-sources:           cpp/hcdt.cpp   extra-libraries:     stdc++   ghc-options:         -Wall -optcxx-std=c++11