diff --git a/cpp/CDT.h b/cpp/CDT.h new file mode 100644 --- /dev/null +++ b/cpp/CDT.h @@ -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 +#include +#include +#include +#include +#include +#include + +/// 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 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 mapping; ///< vertex index mapping + std::vector 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 +void RemoveDuplicates( + std::vector& vertices, + const std::vector& 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 +CDT_EXPORT DuplicatesInfo RemoveDuplicates(std::vector >& 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& 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& edges, const std::vector& 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& 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 +CDT_EXPORT DuplicatesInfo RemoveDuplicatesAndRemapEdges( + std::vector >& vertices, + std::vector& 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 +EdgeToPiecesMapping(const unordered_map& 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 +CDT_EXPORT unordered_map > EdgeToSplitVertices( + const unordered_map& edgeToPieces, + const std::vector >& vertices); + +/// @} + +/// @} + +} // namespace CDT + +//***************************************************************************** +// Implementations of template functionlity +//***************************************************************************** +// hash for CDT::V2d +#ifdef CDT_CXX11_IS_SUPPORTED +namespace std +#else +namespace boost +#endif +{ +template +struct hash > +{ + size_t operator()(const CDT::V2d& xy) const + { +#ifdef CDT_CXX11_IS_SUPPORTED + typedef std::hash Hasher; +#else + typedef boost::hash 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, std::size_t> PosToIndex; + PosToIndex uniqueVerts; + const std::size_t verticesSize = std::distance(first, last); + DuplicatesInfo di = { + std::vector(verticesSize), std::vector()}; + 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::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 +void RemoveDuplicates( + std::vector& vertices, + const std::vector& 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& mapping, + TGetEdgeVertexStart getStart, + TGetEdgeVertexEnd getEnd, + TMakeEdgeFromStartAndEnd makeEdge) +{ + for(; first != last; ++first) + { + *first = makeEdge( + static_cast(mapping[getStart(*first)]), + static_cast(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& vertices, + TGetVertexCoordX getX, + TGetVertexCoordY getY, + const TEdgeIter edgesFirst, + const TEdgeIter edgesLast, + TGetEdgeVertexStart getStart, + TGetEdgeVertexEnd getEnd, + TMakeEdgeFromStartAndEnd makeEdge) +{ + const DuplicatesInfo di = + FindDuplicates(vertices.begin(), vertices.end(), getX, getY); + RemoveDuplicates(vertices, di.duplicates); + RemapEdges(edgesFirst, edgesLast, di.mapping, getStart, getEnd, makeEdge); + return di; +} + +template +unordered_map > EdgeToSplitVertices( + const unordered_map& edgeToPieces, + const std::vector >& vertices) +{ + typedef std::pair VertCoordPair; + struct ComparePred + { + bool operator()(const VertCoordPair& a, const VertCoordPair& b) const + { + return a.second < b.second; + } + } comparePred; + + unordered_map > edgeToSplitVerts; + typedef unordered_map::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 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 vv = {it->v1(), it->v2()}; + typedef array::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 > val = + std::make_pair(e, std::vector()); + val.second.reserve(splitVerts.size()); + typedef typename std::vector::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 diff --git a/cpp/CDT.hpp b/cpp/CDT.hpp new file mode 100644 --- /dev/null +++ b/cpp/CDT.hpp @@ -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 +#include +#include +#include + +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 +DuplicatesInfo RemoveDuplicates(std::vector >& vertices) +{ + const DuplicatesInfo di = FindDuplicates( + vertices.begin(), vertices.end(), getX_V2d, getY_V2d); + RemoveDuplicates(vertices, di.duplicates); + return di; +} + +CDT_INLINE_IF_HEADER_ONLY void +RemapEdges(std::vector& edges, const std::vector& mapping) +{ + RemapEdges( + edges.begin(), + edges.end(), + mapping, + edge_get_v1, + edge_get_v2, + edge_make); +} + +template +DuplicatesInfo RemoveDuplicatesAndRemapEdges( + std::vector >& vertices, + std::vector& edges) +{ + return RemoveDuplicatesAndRemapEdges( + vertices, + getX_V2d, + getY_V2d, + 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 +EdgeToPiecesMapping(const unordered_map& pieceToOriginals) +{ + unordered_map originalToPieces; + typedef unordered_map::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 diff --git a/cpp/CDTUtils.h b/cpp/CDTUtils.h new file mode 100644 --- /dev/null +++ b/cpp/CDTUtils.h @@ -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 +#include +#include +#include + +#ifdef CDT_USE_BOOST +#include +#endif + +// use fall-backs for c++11 features +#ifdef CDT_CXX11_IS_SUPPORTED + +#include +#include +#include +#include +#include +#include +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 +#include +#include +#include +#include +#include +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 +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 +const T& getX_V2d(const V2d& v) +{ + return v.x; +} + +/// Y-coordinate getter for V2d +template +const T& getY_V2d(const V2d& v) +{ + return v.y; +} + +/// If two 2D vectors are exactly equal +template +bool operator==(const CDT::V2d& lhs, const CDT::V2d& 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::max()); +/// Constant representing no valid vertex for a triangle +const static VertInd noVertex(std::numeric_limits::max()); + +typedef std::vector TriIndVec; ///< Vector of triangle indices +typedef array VerticesArr3; ///< array of three vertex indices +typedef array NeighborsArr3; ///< array of three neighbors + +/// 2D bounding box +template +struct CDT_EXPORT Box2d +{ + V2d min; ///< min box corner + V2d max; ///< max box corner + + /// Envelop box around a point + void envelopPoint(const V2d& 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 envelopBox( + TVertexIter first, + TVertexIter last, + TGetVertexCoordX getX, + TGetVertexCoordY getY) +{ + const T max = std::numeric_limits::max(); + Box2d 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 +CDT_EXPORT Box2d envelopBox(const std::vector >& 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& verts() const; + +private: + std::pair 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 EdgeVec; ///< Vector of edges +typedef unordered_set EdgeUSet; ///< Hash table of edges +typedef unordered_set TriIndUSet; ///< Hash table of triangles +typedef unordered_map TriIndUMap; ///< Triangle hash map +#ifdef CDT_USE_BOOST +/// Flat hash table of triangles +typedef boost::container::flat_set 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& vertices, const array& neighbors) + { + Triangle t = {vertices, neighbors}; + return t; + } +}; + +typedef std::vector 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 +CDT_EXPORT T orient2D(const V2d& p, const V2d& v1, const V2d& v2); + +/// Check if point lies to the left of, to the right of, or on a line +template +CDT_EXPORT PtLineLocation::Enum locatePointLine( + const V2d& p, + const V2d& v1, + const V2d& v2, + T orientationTolerance = T(0)); + +/// Classify value of orient2d predicate +template +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 +CDT_EXPORT PtTriLocation::Enum locatePointTriangle( + const V2d& p, + const V2d& v1, + const V2d& v2, + const V2d& 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 +CDT_EXPORT bool isInCircumcircle( + const V2d& p, + const V2d& v1, + const V2d& v2, + const V2d& 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 +CDT_EXPORT T distance(const V2d& a, const V2d& b); + +/// Squared distance between two 2D points +template +CDT_EXPORT T distanceSquared(const V2d& a, const V2d& 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 +{ + /// Hash operator + std::size_t operator()(const CDT::VertInd& vi) const + { + return std::hash()(vi.t); + } +}; + +/// Triangle index hasher +template <> +struct hash +{ + /// Hash operator + std::size_t operator()(const CDT::TriInd& vi) const + { + return std::hash()(vi.t); + } +}; + +#endif // CDT_USE_STRONG_TYPING + +/// Edge hasher +template <> +struct hash +{ + /// 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 Hasher; +#else + typedef boost::hash Hasher; +#endif + seed ^= Hasher()(key) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } + static std::size_t hashEdge(const CDT::Edge& e) + { + const std::pair& 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 diff --git a/cpp/CDTUtils.hpp b/cpp/CDTUtils.hpp new file mode 100644 --- /dev/null +++ b/cpp/CDTUtils.hpp @@ -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 + +namespace CDT +{ + +//***************************************************************************** +// V2d +//***************************************************************************** +template +V2d V2d::make(const T x, const T y) +{ + V2d out = {x, y}; + return out; +} + +//***************************************************************************** +// Box2d +//***************************************************************************** +template +Box2d envelopBox(const std::vector >& vertices) +{ + return envelopBox( + vertices.begin(), vertices.end(), getX_V2d, getY_V2d); +} + +//***************************************************************************** +// 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& 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(location - PtTriLocation::OnEdge1); +} + +template +T orient2D(const V2d& p, const V2d& v1, const V2d& v2) +{ + return predicates::adaptive::orient2d(v1.x, v1.y, v2.x, v2.y, p.x, p.y); +} + +template +PtLineLocation::Enum locatePointLine( + const V2d& p, + const V2d& v1, + const V2d& v2, + const T orientationTolerance) +{ + return classifyOrientation(orient2D(p, v1, v2), orientationTolerance); +} + +template +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 +PtTriLocation::Enum locatePointTriangle( + const V2d& p, + const V2d& v1, + const V2d& v2, + const V2d& 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 +bool isInCircumcircle( + const V2d& p, + const V2d& v1, + const V2d& v2, + const V2d& 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 +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 +T distance(const T ax, const T ay, const T bx, const T by) +{ + return std::sqrt(distanceSquared(ax, ay, bx, by)); +} + +template +T distance(const V2d& a, const V2d& b) +{ + return distance(a.x, a.y, b.x, b.y); +} + +template +T distanceSquared(const V2d& a, const V2d& b) +{ + return distanceSquared(a.x, a.y, b.x, b.y); +} + +} // namespace CDT diff --git a/cpp/KDTree.h b/cpp/KDTree.h new file mode 100644 --- /dev/null +++ b/cpp/KDTree.h @@ -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 + +#ifndef KDTREE_KDTREE_H +#define KDTREE_KDTREE_H + +#include "CDTUtils.h" + +#include +#include + +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 point_type; + typedef CDT::VertInd point_index; + typedef std::pair value_type; + typedef std::vector point_data_vec; + typedef point_data_vec::const_iterator pd_cit; + typedef CDT::VertInd node_index; + typedef CDT::array 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::max(), + -std::numeric_limits::max())) + , m_max(point_type::make( + std::numeric_limits::max(), + std::numeric_limits::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& 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& points) const + { + value_type out; + int iTask = -1; + coord_type minDistSq = std::numeric_limits::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(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(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( + 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& 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 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 m_tasksStack; +}; + +} // namespace KDTree + +#endif // header guard diff --git a/cpp/LocatorKDTree.h b/cpp/LocatorKDTree.h new file mode 100644 --- /dev/null +++ b/cpp/LocatorKDTree.h @@ -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 >& points) + { + typedef V2d V2d_t; + V2d_t min = points.front(); + V2d_t max = min; + typedef typename std::vector::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 >& points) + { + m_kdTree.insert(i, points); + } + /// Find nearest point using R-tree + VertInd nearPoint( + const V2d& pos, + const std::vector >& 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 diff --git a/cpp/Triangulation.h b/cpp/Triangulation.h new file mode 100644 --- /dev/null +++ b/cpp/Triangulation.h @@ -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 +#include +#include +#include +#include +#include +#include + +/// 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 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 > +class CDT_EXPORT Triangulation +{ +public: + typedef std::vector > 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 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 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 >& 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. + * No new vertices are inserted. + * @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. + * Make sure there are no erroneous duplicates. + * @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. + * No new vertices are inserted. + * @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. + * Make sure there are no erroneous duplicates. + * @tparam edges constraint edges + */ + void insertEdges(const std::vector& 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. New vertices are inserted. + * @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. + * Make sure there are no erroneous duplicates. + * @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. New vertices are inserted. + * @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. + * Make sure there are no erroneous duplicates. + * @tparam edges edges to conform to + */ + void conformToEdges(const std::vector& 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 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& box); + void addNewVertex(const V2d& pos, const TriIndVec& tris); + void insertVertex(VertInd iVert); + void ensureDelaunayByEdgeFlips( + const V2d& v, + VertInd iVert, + std::stack& triStack); + /// Flip fixed edges and return a list of flipped fixed edges + std::vector 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 intersectedTriangle( + VertInd iA, + const std::vector& candidates, + const V2d& a, + const V2d& b, + T orientationTolerance = T(0)) const; + /// Returns indices of three resulting triangles + std::stack insertPointInTriangle(VertInd v, TriInd iT); + /// Returns indices of four resulting triangles + std::stack insertPointOnEdge(VertInd v, TriInd iT1, TriInd iT2); + array trianglesAt(const V2d& pos) const; + array walkingSearchTrianglesAt(const V2d& pos) const; + TriInd walkTriangles(VertInd startVertex, const V2d& pos) const; + bool isFlipNeeded( + const V2d& v, + VertInd iV, + VertInd iV1, + VertInd iV2, + VertInd iV3) const; + bool + isFlipNeeded(const V2d& 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::const_iterator pointsFirst, + std::vector::const_iterator pointsLast); + VertInd findDelaunayPoint( + VertInd ia, + VertInd ib, + std::vector::const_iterator pointsFirst, + std::vector::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 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 peelLayer( + std::stack seeds, + LayerDepth layerDepth, + std::vector& triDepths) const; + + std::vector 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 +void random_shuffle(RandomIt first, RandomIt last) +{ + typename std::iterator_traits::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 +template < + typename TVertexIter, + typename TGetVertexCoordX, + typename TGetVertexCoordY> +void Triangulation::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(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::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 ii(std::distance(first, last)); + typedef std::vector::iterator Iter; + VertInd value = static_cast(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 +template < + typename TEdgeIter, + typename TGetEdgeVertexStart, + typename TGetEdgeVertexEnd> +void Triangulation::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 +template < + typename TEdgeIter, + typename TGetEdgeVertexStart, + typename TGetEdgeVertexEnd> +void Triangulation::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 diff --git a/cpp/Triangulation.hpp b/cpp/Triangulation.hpp new file mode 100644 --- /dev/null +++ b/cpp/Triangulation.hpp @@ -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 +#include +#include +#include + +namespace CDT +{ + +typedef std::deque TriDeque; + +namespace detail +{ + +/// Needed for c++03 compatibility (no uniform initialization available) +template +array arr3(const T& v0, const T& v1, const T& v2) +{ + const array 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 +Triangulation::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 +Triangulation::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 +Triangulation::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 +Triangulation::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 +void Triangulation::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 +void Triangulation::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(); +} + +template +void Triangulation::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 +void Triangulation::eraseOuterTriangles() +{ + // make dummy triangles adjacent to super-triangle's vertices + const std::stack seed(std::deque(1, vertTris[0].front())); + const TriIndUSet toErase = growToBoundary(seed); + finalizeTriangulation(toErase); +} + +template +void Triangulation::eraseOuterTrianglesAndHoles() +{ + const std::vector 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(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 +void Triangulation::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 +void Triangulation::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 updatedOverlapCount; + typedef unordered_map::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 updatedPieceToOriginals; + typedef unordered_map::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 +void Triangulation::initializedWithCustomSuperGeometry() +{ + m_nearPtLocator.initialize(vertices); + m_nTargetVerts = vertices.size(); + m_superGeomType = SuperGeometryType::Custom; +} + +template +TriIndUSet Triangulation::growToBoundary( + std::stack 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 +void Triangulation::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 +TriInd Triangulation::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 +TriInd Triangulation::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 +void Triangulation::insertEdges( + const std::vector& edges) +{ + insertEdges(edges.begin(), edges.end(), edge_get_v1, edge_get_v2); +} + +template +void Triangulation::conformToEdges( + const std::vector& edges) +{ + conformToEdges(edges.begin(), edges.end(), edge_get_v1, edge_get_v2); +} + +template +void Triangulation::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 +void insert_unique(std::vector& 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 +void insert_unique( + std::vector& to, + const std::vector& from) +{ + typedef typename std::vector::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 +void Triangulation::fixEdge( + const Edge& edge, + const Edge& originalEdge) +{ + fixEdge(edge); + if(edge != originalEdge) + detail::insert_unique(pieceToOriginals[edge], originalEdge); +} + +template +void Triangulation::fixEdge( + const Edge& edge, + const BoundaryOverlapCount overlaps) +{ + fixedEdges.insert(edge); + overlapCount[edge] = overlaps; // override overlap counter +} + +namespace detail +{ + +template +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 +V2d intersectionPosition( + const V2d& a, + const V2d& b, + const V2d& c, + const V2d& 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::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::make(lerp(c.x, d.x, t), lerp(c.y, d.y, t)); + } +} + +} // namespace detail + +template +void Triangulation::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& a = vertices[iA]; + const V2d& 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 intersected(1, iT); + std::vector ptsLeft(1, iVleft); + std::vector 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 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(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::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 newV = detail::intersectionPosition( + vertices[iA], + vertices[iB], + vertices[iVleft], + vertices[iVright]); + addNewVertex(newV, TriIndVec()); + std::stack 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::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 +void Triangulation::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& a = vertices[iA]; + const V2d& 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 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(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::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 newV = detail::intersectionPosition( + vertices[iA], + vertices[iB], + vertices[iVleft], + vertices[iVright]); + addNewVertex(newV, TriIndVec()); + std::stack 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(vertices.size()); + const V2d& start = vertices[iA]; + const V2d& end = vertices[iB]; + addNewVertex( + V2d::make((start.x + end.x) / T(2), (start.y + end.y) / T(2)), + TriIndVec()); + const std::vector 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::const_iterator it = flippedFixedEdges.begin(); + it != flippedFixedEdges.end(); + ++it) + { + fixedEdges.erase(*it); + + BoundaryOverlapCount prevOverlaps = 0; + const unordered_map::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::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 +tuple +Triangulation::intersectedTriangle( + const VertInd iA, + const std::vector& candidates, + const V2d& a, + const V2d& b, + const T orientationTolerance) const +{ + typedef std::vector::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 +void Triangulation::addSuperTriangle(const Box2d& box) +{ + m_nTargetVerts = 3; + m_superGeomType = SuperGeometryType::SuperTriangle; + + const V2d 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 posV1 = {center.x - shiftX, center.y - r}; + const V2d posV2 = {center.x + shiftX, center.y - r}; + const V2d 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 +void Triangulation::addNewVertex( + const V2d& pos, + const TriIndVec& tris) +{ + vertices.push_back(pos); + vertTris.push_back(tris); +} + +template +std::vector +Triangulation::insertVertex_FlipFixedEdges( + const VertInd iVert) +{ + std::vector flippedFixedEdges; + + const V2d& v = vertices[iVert]; + array trisAt = walkingSearchTrianglesAt(v); + std::stack 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 +void Triangulation::insertVertex(const VertInd iVert) +{ + const V2d& v = vertices[iVert]; + array trisAt = walkingSearchTrianglesAt(v); + std::stack triStack = + trisAt[1] == noNeighbor + ? insertPointInTriangle(iVert, trisAt[0]) + : insertPointOnEdge(iVert, trisAt[0], trisAt[1]); + ensureDelaunayByEdgeFlips(v, iVert, triStack); + m_nearPtLocator.addPoint(iVert, vertices); +} + +template +void Triangulation::ensureDelaunayByEdgeFlips( + const V2d& v, + const VertInd iVert, + std::stack& 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 +bool Triangulation::isFlipNeeded( + const V2d& v, + const VertInd iV, + const VertInd iV1, + const VertInd iV2, + const VertInd iV3) const +{ + const V2d& v1 = vertices[iV1]; + const V2d& v2 = vertices[iV2]; + const V2d& 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 +bool Triangulation::isFlipNeeded( + const V2d& 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 +std::stack Triangulation::insertPointInTriangle( + const VertInd v, + const TriInd iT) +{ + const TriInd iNewT1 = addTriangle(); + const TriInd iNewT2 = addTriangle(); + + Triangle& t = triangles[iT]; + const array vv = t.vertices; + const array 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 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 +std::stack Triangulation::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 newTriangles; + newTriangles.push(iT1); + newTriangles.push(iTnew2); + newTriangles.push(iT2); + newTriangles.push(iTnew1); + return newTriangles; +} + +template +array +Triangulation::trianglesAt(const V2d& pos) const +{ + array out = {noNeighbor, noNeighbor}; + for(TriInd i = TriInd(0); i < TriInd(triangles.size()); ++i) + { + const Triangle& t = triangles[i]; + const V2d& v1 = vertices[t.vertices[0]]; + const V2d& v2 = vertices[t.vertices[1]]; + const V2d& 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 +TriInd Triangulation::walkTriangles( + const VertInd startVertex, + const V2d& 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& vStart = vertices[t.vertices[i]]; + const V2d& 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 +array Triangulation::walkingSearchTrianglesAt( + const V2d& pos) const +{ + array 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& v1 = vertices[t.vertices[0]]; + const V2d& v2 = vertices[t.vertices[1]]; + const V2d& 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 +void Triangulation::flipEdge( + const TriInd iT, + const TriInd iTopo) +{ + Triangle& t = triangles[iT]; + Triangle& tOpo = triangles[iTopo]; + const array& triNs = t.neighbors; + const array& triOpoNs = tOpo.neighbors; + const array& triVs = t.vertices; + const array& 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 +void Triangulation::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 +void Triangulation::addAdjacentTriangle( + const VertInd iVertex, + const TriInd iTriangle) +{ + vertTris[iVertex].push_back(iTriangle); +} + +template +void Triangulation::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 +void Triangulation::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 +void Triangulation::removeAdjacentTriangle( + const VertInd iVertex, + const TriInd iTriangle) +{ + std::vector& tris = vertTris[iVertex]; + tris.erase(std::find(tris.begin(), tris.end(), iTriangle)); +} + +template +TriInd Triangulation::triangulatePseudopolygon( + const VertInd ia, + const VertInd ib, + const std::vector::const_iterator pointsFirst, + const std::vector::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::const_iterator newLast = pointsFirst; + while(*newLast != ic) + ++newLast; + const std::vector::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 +VertInd Triangulation::findDelaunayPoint( + const VertInd ia, + const VertInd ib, + const std::vector::const_iterator pointsFirst, + const std::vector::const_iterator pointsLast) const +{ + assert(pointsFirst != pointsLast); + const V2d& a = vertices[ia]; + const V2d& b = vertices[ib]; + VertInd ic = *pointsFirst; + V2d c = vertices[ic]; + typedef std::vector::const_iterator CIt; + for(CIt it = pointsFirst + 1; it != pointsLast; ++it) + { + const V2d v = vertices[*it]; + if(!isInCircumcircle(v, a, b, c)) + continue; + ic = *it; + c = vertices[ic]; + } + return ic; +} + +template +TriInd Triangulation::pseudopolyOuterTriangle( + const VertInd ia, + const VertInd ib) const +{ + const std::vector& aTris = vertTris[ia]; + const std::vector& bTris = vertTris[ib]; + typedef std::vector::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 +void Triangulation::insertVertices( + const std::vector >& newVertices) +{ + return insertVertices( + newVertices.begin(), newVertices.end(), getX_V2d, getY_V2d); +} + +template +bool Triangulation::isFinalized() const +{ + return vertTris.empty() && !vertices.empty(); +} + +template +unordered_map +Triangulation::peelLayer( + std::stack seeds, + const LayerDepth layerDepth, + std::vector& triDepths) const +{ + unordered_map 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::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 +std::vector +Triangulation::calculateTriangleDepths() const +{ + std::vector triDepths( + triangles.size(), std::numeric_limits::max()); + std::stack seeds(TriDeque(1, vertTris[0].front())); + LayerDepth layerDepth = 0; + LayerDepth deepestSeedDepth = 0; + + unordered_map seedsByDepth; + do + { + const unordered_map& newSeeds = + peelLayer(seeds, layerDepth, triDepths); + + seedsByDepth.erase(layerDepth); + typedef unordered_map::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( + TriDeque(nextLayerSeeds.begin(), nextLayerSeeds.end())); + ++layerDepth; + } while(!seeds.empty() || deepestSeedDepth > layerDepth); + + return triDepths; +} + +} // namespace CDT diff --git a/cpp/hcdt.hpp b/cpp/hcdt.hpp new file mode 100644 --- /dev/null +++ b/cpp/hcdt.hpp @@ -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; diff --git a/cpp/predicates.h b/cpp/predicates.h new file mode 100644 --- /dev/null +++ b/cpp/predicates.h @@ -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 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 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 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 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 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 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 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 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 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 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 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 T insphere(T const*const pa, T const*const pb, T const*const pc, T const*const pd, T const*const pe); + } +} + +#include //abs, fma +#include +#include //pair +#include //accumulate +#include //transform, copy_n, merge +#include //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 +#include // 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 + 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 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 + class Expansion : private ExpansionBase, public stdx::array { + private: + public: + size_t m_size; + template friend class ExpansionBase;//access for base class + template friend class Expansion;//access for expansions of different size + + Expansion() : m_size(0) {} + template Expansion& operator=(const Expansion& 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::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::operator[](m_size++) = v;} + + public: + //estimates of expansion value + T mostSignificant() const {return empty() ? T(0) : stdx::array::operator[](m_size - 1);} + T estimate() const {return std::accumulate(stdx::array::cbegin(), stdx::array::cbegin() + size(), T(0));} + + template Expansion operator+(const Expansion& f) const { + Expansion h; + h.m_size = ExpansionBase::ExpansionSum(this->data(), this->size(), f.data(), f.size(), h.data()); + return h; + } + + void negate() {std::transform(stdx::array::cbegin(), stdx::array::cbegin() + size(), stdx::array::begin(), std::negate());} + Expansion operator-() const {Expansion e = *this; e.negate(); return e;} + template Expansion operator-(const Expansion& f) const {return operator+(-f);} + + Expansion operator*(const T b) const { + Expansion h; + h.m_size = ExpansionBase::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 struct use_fma {static const bool value = (std::is_same::value && fp_fast_fmaf) || + (std::is_same::value && fp_fast_fma) || + (std::is_same::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 bool absLess(const T& a, const T& b) {return std::abs(a) < std::abs(b);} + + template + class ExpansionBase { + private: + static const T Splitter; + + PREDICATES_PORTABLE_STATIC_ASSERT(std::numeric_limits::is_iec559, Requires_IEC_559_IEEE_754_floating_point_type); + PREDICATES_PORTABLE_STATIC_ASSERT(2 == std::numeric_limits::radix, Requires_base_2_floating_point_type); + + //combine result + roundoff error into expansion + static inline Expansion MakeExpansion(const T value, const T tail) { + Expansion 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); + 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 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 Split(const T a) { + const T c = a * Splitter; + const T aBig = c - a; + const T aHi = c - aBig; + return std::pair(aHi, a - aHi); + } + + //roundoff error of x = a * b via dekkers product + static inline T DekkersProduct(const T /*a*/, const std::pair aSplit, const T /*b*/, const std::pair 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 static typename std::enable_if< use_fma::value, S>::type MultTail(const T a, const T b, const T p) {return std::fma(a, b, -p);} + template static typename std::enable_if::value, S>::type MultTail(const T a, const T b, const T p) {return DekkersProduct(a, Split(a), b, Split(b), p);} + + template static typename std::enable_if< use_fma::value, S>::type MultTailPreSplit(const T a, const T b, const std::pair /*bSplit*/, const T p) {return std::fma(a, b, -p);} + template static typename std::enable_if::value, S>::type MultTailPreSplit(const T a, const T b, const std::pair 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 bSplit, const T p) {return DekkersProduct(a, Split(a), b, bSplit, p);} +#endif + //expand a + b + static inline Expansion 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 Minus(const T a, const T b) {return Plus(a, -b);} + + //expand a * b + static inline Expansion 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 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 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 TwoTwoDiffZeroCheck(const T ax, const T by, const T ay, const T bx) { + Expansion 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 ThreeProd(const T a, const T b, const T c) {return (T(0) == a || T(0) == b || T(0) == c) ? Expansion() : Mult(a, b) * c;} + }; + + template const T ExpansionBase::Splitter = static_cast( +#ifdef PREDICATES_CXX11_IS_SUPPORTED + std::exp2((std::numeric_limits::digits + std::numeric_limits::digits%2)/2 + 1) +#else + std::ldexp(T(1), (std::numeric_limits::digits + std::numeric_limits::digits%2)/2 + 1) +#endif + ); +} + + namespace exact { + template T orient2d(T const ax, T const ay, T const bx, T const by, T const cx, T const cy) + { + const detail::Expansion aterms = detail::ExpansionBase::TwoTwoDiff(ax, by, ax, cy); + const detail::Expansion bterms = detail::ExpansionBase::TwoTwoDiff(bx, cy, bx, ay); + const detail::Expansion cterms = detail::ExpansionBase::TwoTwoDiff(cx, ay, cx, by); + const detail::Expansion w = aterms + bterms + cterms; + return w.mostSignificant(); + } + + template 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 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 ab = detail::ExpansionBase::TwoTwoDiff(ax, by, bx, ay); + const detail::Expansion bc = detail::ExpansionBase::TwoTwoDiff(bx, cy, cx, by); + const detail::Expansion cd = detail::ExpansionBase::TwoTwoDiff(cx, dy, dx, cy); + const detail::Expansion da = detail::ExpansionBase::TwoTwoDiff(dx, ay, ax, dy); + const detail::Expansion ac = detail::ExpansionBase::TwoTwoDiff(ax, cy, cx, ay); + const detail::Expansion bd = detail::ExpansionBase::TwoTwoDiff(bx, dy, dx, by); + + const detail::Expansion abc = ab + bc - ac; + const detail::Expansion bcd = bc + cd - bd; + const detail::Expansion cda = cd + da + ac; + const detail::Expansion dab = da + ab + bd; + + const detail::Expansion adet = bcd * ax * ax + bcd * ay * ay; + const detail::Expansion bdet = cda * bx * -bx + cda * by * -by; + const detail::Expansion cdet = dab * cx * cx + dab * cy * cy; + const detail::Expansion ddet = abc * dx * -dx + abc * dy * -dy; + + const detail::Expansion deter = (adet + bdet) + (cdet + ddet); + return deter.mostSignificant(); + } + + template 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 T orient3d(T const*const pa, T const*const pb, T const*const pc, T const*const pd) { + const detail::Expansion ab = detail::ExpansionBase::TwoTwoDiff(pa[0], pb[1], pb[0], pa[1]); + const detail::Expansion bc = detail::ExpansionBase::TwoTwoDiff(pb[0], pc[1], pc[0], pb[1]); + const detail::Expansion cd = detail::ExpansionBase::TwoTwoDiff(pc[0], pd[1], pd[0], pc[1]); + const detail::Expansion da = detail::ExpansionBase::TwoTwoDiff(pd[0], pa[1], pa[0], pd[1]); + const detail::Expansion ac = detail::ExpansionBase::TwoTwoDiff(pa[0], pc[1], pc[0], pa[1]); + const detail::Expansion bd = detail::ExpansionBase::TwoTwoDiff(pb[0], pd[1], pd[0], pb[1]); + + const detail::Expansion abc = ab + bc - ac; + const detail::Expansion bcd = bc + cd - bd; + const detail::Expansion cda = cd + da + ac; + const detail::Expansion dab = da + ab + bd; + + const detail::Expansion adet = bcd * pa[2]; + const detail::Expansion bdet = cda * -pb[2]; + const detail::Expansion cdet = dab * pc[2]; + const detail::Expansion ddet = abc * -pd[2]; + + const detail::Expansion 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 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 ab = detail::ExpansionBase::TwoTwoDiff(pa[0], pb[1], pb[0], pa[1]); + const detail::Expansion bc = detail::ExpansionBase::TwoTwoDiff(pb[0], pc[1], pc[0], pb[1]); + const detail::Expansion cd = detail::ExpansionBase::TwoTwoDiff(pc[0], pd[1], pd[0], pc[1]); + const detail::Expansion de = detail::ExpansionBase::TwoTwoDiff(pd[0], pe[1], pe[0], pd[1]); + const detail::Expansion ea = detail::ExpansionBase::TwoTwoDiff(pe[0], pa[1], pa[0], pe[1]); + const detail::Expansion ac = detail::ExpansionBase::TwoTwoDiff(pa[0], pc[1], pc[0], pa[1]); + const detail::Expansion bd = detail::ExpansionBase::TwoTwoDiff(pb[0], pd[1], pd[0], pb[1]); + const detail::Expansion ce = detail::ExpansionBase::TwoTwoDiff(pc[0], pe[1], pe[0], pc[1]); + const detail::Expansion da = detail::ExpansionBase::TwoTwoDiff(pd[0], pa[1], pa[0], pd[1]); + const detail::Expansion eb = detail::ExpansionBase::TwoTwoDiff(pe[0], pb[1], pb[0], pe[1]); + + const detail::Expansion abc = bc * pa[2] + ac * -pb[2] + ab * pc[2]; + const detail::Expansion bcd = cd * pb[2] + bd * -pc[2] + bc * pd[2]; + const detail::Expansion cde = de * pc[2] + ce * -pd[2] + cd * pe[2]; + const detail::Expansion dea = ea * pd[2] + da * -pe[2] + de * pa[2]; + const detail::Expansion eab = ab * pe[2] + eb * -pa[2] + ea * pb[2]; + const detail::Expansion abd = bd * pa[2] + da * pb[2] + ab * pd[2]; + const detail::Expansion bce = ce * pb[2] + eb * pc[2] + bc * pe[2]; + const detail::Expansion cda = da * pc[2] + ac * pd[2] + cd * pa[2]; + const detail::Expansion deb = eb * pd[2] + bd * pe[2] + de * pb[2]; + const detail::Expansion eac = ac * pe[2] + ce * pa[2] + ea * pc[2]; + + const detail::Expansion bcde = (cde + bce) - (deb + bcd); + const detail::Expansion cdea = (dea + cda) - (eac + cde); + const detail::Expansion deab = (eab + deb) - (abd + dea); + const detail::Expansion eabc = (abc + eac) - (bce + eab); + const detail::Expansion abcd = (bcd + abd) - (cda + abc); + + const detail::Expansion adet = bcde * pa[0] * pa[0] + bcde * pa[1] * pa[1] + bcde * pa[2] * pa[2]; + const detail::Expansion bdet = cdea * pb[0] * pb[0] + cdea * pb[1] * pb[1] + cdea * pb[2] * pb[2]; + const detail::Expansion cdet = deab * pc[0] * pc[0] + deab * pc[1] * pc[1] + deab * pc[2] * pc[2]; + const detail::Expansion ddet = eabc * pd[0] * pd[0] + eabc * pd[1] * pd[1] + eabc * pd[2] * pd[2]; + const detail::Expansion edet = abcd * pe[0] * pe[0] + abcd * pe[1] * pe[1] + abcd * pe[2] * pe[2]; + + const detail::Expansion deter = (adet + bdet) + ((cdet + ddet) + edet); + return deter.mostSignificant(); + } + } + + template + const T& Epsilon() + { + static const T epsilon = static_cast( +#ifdef PREDICATES_CXX11_IS_SUPPORTED + std::exp2(-std::numeric_limits::digits) +#else + std::ldexp(T(1), -std::numeric_limits::digits) +#endif + ); + return epsilon; + } + + template + 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 const T Constants::epsilon = Epsilon(); + template const T Constants::resulterrbound = (T( 3) + T( 8) * Epsilon()) * Epsilon(); + template const T Constants::ccwerrboundA = (T( 3) + T( 16) * Epsilon()) * Epsilon(); + template const T Constants::ccwerrboundB = (T( 2) + T( 12) * Epsilon()) * Epsilon(); + template const T Constants::ccwerrboundC = (T( 9) + T( 64) * Epsilon()) * Epsilon() * Epsilon(); + template const T Constants::o3derrboundA = (T( 7) + T( 56) * Epsilon()) * Epsilon(); + template const T Constants::o3derrboundB = (T( 3) + T( 28) * Epsilon()) * Epsilon(); + template const T Constants::o3derrboundC = (T(26) + T( 288) * Epsilon()) * Epsilon() * Epsilon(); + template const T Constants::iccerrboundA = (T(10) + T( 96) * Epsilon()) * Epsilon(); + template const T Constants::iccerrboundB = (T( 4) + T( 48) * Epsilon()) * Epsilon(); + template const T Constants::iccerrboundC = (T(44) + T( 576) * Epsilon()) * Epsilon() * Epsilon(); + template const T Constants::isperrboundA = (T(16) + T( 224) * Epsilon()) * Epsilon(); + template const T Constants::isperrboundB = (T( 5) + T( 72) * Epsilon()) * Epsilon(); + template const T Constants::isperrboundC = (T(71) + T(1408) * Epsilon()) * Epsilon() * Epsilon(); + + namespace adaptive { + template 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::ccwerrboundA * detsum; + if(std::abs(det) >= std::abs(errbound)) return det; + + const detail::Expansion B = detail::ExpansionBase::TwoTwoDiff(acx, bcy, acy, bcx); + det = B.estimate(); + errbound = Constants::ccwerrboundB * detsum; + if(std::abs(det) >= std::abs(errbound)) return det; + + const T acxtail = detail::ExpansionBase::MinusTail(ax, cx, acx); + const T bcxtail = detail::ExpansionBase::MinusTail(bx, cx, bcx); + const T acytail = detail::ExpansionBase::MinusTail(ay, cy, acy); + const T bcytail = detail::ExpansionBase::MinusTail(by, cy, bcy); + if(T(0) == acxtail && T(0) == bcxtail && T(0) == acytail && T(0) == bcytail) return det; + + errbound = Constants::ccwerrboundC * detsum + Constants::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 D = ((B + detail::ExpansionBase::TwoTwoDiff(acxtail, bcy, acytail, bcx)) + detail::ExpansionBase::TwoTwoDiff(acx, bcytail, acy, bcxtail)) + detail::ExpansionBase::TwoTwoDiff(acxtail, bcytail, acytail, bcxtail); + return D.mostSignificant(); + } + + template 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 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::iccerrboundA * permanent; + if(std::abs(det) >= std::abs(errbound)) return det; + + const detail::Expansion bc = detail::ExpansionBase::TwoTwoDiff(bdx, cdy, cdx, bdy); + const detail::Expansion ca = detail::ExpansionBase::TwoTwoDiff(cdx, ady, adx, cdy); + const detail::Expansion ab = detail::ExpansionBase::TwoTwoDiff(adx, bdy, bdx, ady); + const detail::Expansion adet = bc * adx * adx + bc * ady * ady; + const detail::Expansion bdet = ca * bdx * bdx + ca * bdy * bdy; + const detail::Expansion cdet = ab * cdx * cdx + ab * cdy * cdy; + const detail::Expansion fin1 = adet + bdet + cdet; + det = fin1.estimate(); + errbound = Constants::iccerrboundB * permanent; + if(std::abs(det) >= std::abs(errbound)) return det; + + const T adxtail = detail::ExpansionBase::MinusTail(ax, dx, adx); + const T adytail = detail::ExpansionBase::MinusTail(ay, dy, ady); + const T bdxtail = detail::ExpansionBase::MinusTail(bx, dx, bdx); + const T bdytail = detail::ExpansionBase::MinusTail(by, dy, bdy); + const T cdxtail = detail::ExpansionBase::MinusTail(cx, dx, cdx); + const T cdytail = detail::ExpansionBase::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::iccerrboundC * permanent + Constants::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 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 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::o3derrboundA * permanent; + if(std::abs(det) >= std::abs(errbound)) return det; + + const detail::Expansion bc = detail::ExpansionBase::TwoTwoDiff(bdx, cdy, cdx, bdy); + const detail::Expansion ca = detail::ExpansionBase::TwoTwoDiff(cdx, ady, adx, cdy); + const detail::Expansion ab = detail::ExpansionBase::TwoTwoDiff(adx, bdy, bdx, ady); + const detail::Expansion fin1 = (bc * adz + ca * bdz) + ab * cdz; + det = fin1.estimate(); + errbound = Constants::o3derrboundB * permanent; + if(std::abs(det) >= std::abs(errbound)) return det; + + const T adxtail = detail::ExpansionBase::MinusTail(pa[0], pd[0], adx); + const T bdxtail = detail::ExpansionBase::MinusTail(pb[0], pd[0], bdx); + const T cdxtail = detail::ExpansionBase::MinusTail(pc[0], pd[0], cdx); + const T adytail = detail::ExpansionBase::MinusTail(pa[1], pd[1], ady); + const T bdytail = detail::ExpansionBase::MinusTail(pb[1], pd[1], bdy); + const T cdytail = detail::ExpansionBase::MinusTail(pc[1], pd[1], cdy); + const T adztail = detail::ExpansionBase::MinusTail(pa[2], pd[2], adz); + const T bdztail = detail::ExpansionBase::MinusTail(pb[2], pd[2], bdz); + const T cdztail = detail::ExpansionBase::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::o3derrboundC * permanent + Constants::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 bct = detail::ExpansionBase::TwoTwoDiffZeroCheck(bdxtail, cdy, bdytail, cdx) + detail::ExpansionBase::TwoTwoDiffZeroCheck(cdytail, bdx, cdxtail, bdy); + const detail::Expansion cat = detail::ExpansionBase::TwoTwoDiffZeroCheck(cdxtail, ady, cdytail, adx) + detail::ExpansionBase::TwoTwoDiffZeroCheck(adytail, cdx, adxtail, cdy); + const detail::Expansion abt = detail::ExpansionBase::TwoTwoDiffZeroCheck(adxtail, bdy, adytail, bdx) + detail::ExpansionBase::TwoTwoDiffZeroCheck(bdytail, adx, bdxtail, ady); + const detail::Expansion fin2 = fin1 + bct * adz + cat * bdz + abt * cdz + bc * adztail + ca * bdztail + ab * cdztail + + detail::ExpansionBase::ThreeProd( adxtail, bdytail, cdz) + detail::ExpansionBase::ThreeProd( adxtail, bdytail, cdztail) + + detail::ExpansionBase::ThreeProd(-adxtail, cdytail, bdz) + detail::ExpansionBase::ThreeProd(-adxtail, cdytail, bdztail) + + detail::ExpansionBase::ThreeProd( bdxtail, cdytail, adz) + detail::ExpansionBase::ThreeProd( bdxtail, cdytail, adztail) + + detail::ExpansionBase::ThreeProd(-bdxtail, adytail, cdz) + detail::ExpansionBase::ThreeProd(-bdxtail, adytail, cdztail) + + detail::ExpansionBase::ThreeProd( cdxtail, adytail, bdz) + detail::ExpansionBase::ThreeProd( cdxtail, adytail, bdztail) + + detail::ExpansionBase::ThreeProd(-cdxtail, bdytail, adz) + detail::ExpansionBase::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 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::isperrboundA * permanent; + if(std::abs(det) >= std::abs(errbound)) return det; + } + + const detail::Expansion ab = detail::ExpansionBase::TwoTwoDiff(aex, bey, bex, aey); + const detail::Expansion bc = detail::ExpansionBase::TwoTwoDiff(bex, cey, cex, bey); + const detail::Expansion cd = detail::ExpansionBase::TwoTwoDiff(cex, dey, dex, cey); + const detail::Expansion da = detail::ExpansionBase::TwoTwoDiff(dex, aey, aex, dey); + const detail::Expansion ac = detail::ExpansionBase::TwoTwoDiff(aex, cey, cex, aey); + const detail::Expansion bd = detail::ExpansionBase::TwoTwoDiff(bex, dey, dex, bey); + const detail::Expansion temp24a = bc * dez + (cd * bez + bd * -cez); + const detail::Expansion temp24b = cd * aez + (da * cez + ac * dez); + const detail::Expansion temp24c = da * bez + (ab * dez + bd * aez); + const detail::Expansion temp24d = ab * cez + (bc * aez + ac * -bez); + const detail::Expansion adet = temp24a * aex * -aex + temp24a * aey * -aey + temp24a * aez * -aez; + const detail::Expansion bdet = temp24b * bex * bex + temp24b * bey * bey + temp24b * bez * bez; + const detail::Expansion cdet = temp24c * cex * -cex + temp24c * cey * -cey + temp24c * cez * -cez; + const detail::Expansion ddet = temp24d * dex * dex + temp24d * dey * dey + temp24d * dez * dez; + const detail::Expansion fin1 = (adet + bdet) + (cdet + ddet); + T det = fin1.estimate(); + T errbound = Constants::isperrboundB * permanent; + if(std::abs(det) >= std::abs(errbound)) return det; + + const T aextail = detail::ExpansionBase::MinusTail(pa[0], pe[0], aex); + const T aeytail = detail::ExpansionBase::MinusTail(pa[1], pe[1], aey); + const T aeztail = detail::ExpansionBase::MinusTail(pa[2], pe[2], aez); + const T bextail = detail::ExpansionBase::MinusTail(pb[0], pe[0], bex); + const T beytail = detail::ExpansionBase::MinusTail(pb[1], pe[1], bey); + const T beztail = detail::ExpansionBase::MinusTail(pb[2], pe[2], bez); + const T cextail = detail::ExpansionBase::MinusTail(pc[0], pe[0], cex); + const T ceytail = detail::ExpansionBase::MinusTail(pc[1], pe[1], cey); + const T ceztail = detail::ExpansionBase::MinusTail(pc[2], pe[2], cez); + const T dextail = detail::ExpansionBase::MinusTail(pd[0], pe[0], dex); + const T deytail = detail::ExpansionBase::MinusTail(pd[1], pe[1], dey); + const T deztail = detail::ExpansionBase::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::isperrboundC * permanent + Constants::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 diff --git a/cpp/remove_at.hpp b/cpp/remove_at.hpp new file mode 100644 --- /dev/null +++ b/cpp/remove_at.hpp @@ -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 +#include + +/*! + * Remove elements in the range [first; last) with indices from the sorted + * unique range [ii_first, ii_last) + */ +template +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::difference_type diff_t; + typedef typename std::iterator_traits::value_type ind_t; + ForwardIt destination = first + static_cast(*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(*(ii_first - 1)) + 1; + const ForwardIt source_last = + ii_first != ii_last ? first + static_cast(*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 diff --git a/hcdt.cabal b/hcdt.cabal --- a/hcdt.cabal +++ b/hcdt.cabal @@ -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