From c53e3613042932733b526d97e66252bf8424f371 Mon Sep 17 00:00:00 2001 From: artem-ogre Date: Wed, 20 Oct 2021 22:02:05 +0200 Subject: [PATCH] #40 Performance improvements - Implement randomized vertex insertion option and make it default - Allow custom nearest point locators - Implement boost::rtree and random-nearest as locators - Implement kd-tree for nearest point location - CDT code refactoring, some APIs changed - Make necessary adjustments to visualizer code --- CDT/CMakeLists.txt | 11 +- CDT/extras/VerifyTopology.h | 10 +- CDT/include/CDT.h | 148 +++++++---- CDT/include/CDT.hpp | 389 +++++++++++++++-------------- CDT/include/CDTUtils.h | 48 ++-- CDT/include/CDTUtils.hpp | 76 ++---- CDT/include/KDTree.h | 339 +++++++++++++++++++++++++ CDT/include/LocatorBoostRTree.h | 55 ++++ CDT/include/LocatorKDTree.h | 42 ++++ CDT/include/LocatorNearestRandom.h | 69 +++++ CDT/include/PointRTree.h | 122 --------- CDT/src/CDT.cpp | 26 +- visualizer/main.cpp | 53 ++-- 13 files changed, 896 insertions(+), 492 deletions(-) create mode 100644 CDT/include/KDTree.h create mode 100644 CDT/include/LocatorBoostRTree.h create mode 100644 CDT/include/LocatorKDTree.h create mode 100644 CDT/include/LocatorNearestRandom.h delete mode 100644 CDT/include/PointRTree.h diff --git a/CDT/CMakeLists.txt b/CDT/CMakeLists.txt index b15fd639..a6ef7fa3 100644 --- a/CDT/CMakeLists.txt +++ b/CDT/CMakeLists.txt @@ -26,6 +26,10 @@ option(CDT_USE_BOOST option(CDT_USE_AS_COMPILED_LIBRARY "If enabled templates for float and double will be instantiated and compiled into a library") +option(CDT_USE_64_BIT_INDEX_TYPE + "If enabled 64bits are used to store vertex/triangle index types. Otherwise 32bits are used (up to 4.2bn items)" + OFF) + # check if Boost is needed if(cxx_std_11 IN_LIST CMAKE_CXX_COMPILE_FEATURES) # Work-around as AppleClang 11 defaults to c++98 by default @@ -40,6 +44,7 @@ endif() message(STATUS "CDT_USE_BOOST is ${CDT_USE_BOOST}") message(STATUS "CDT_USE_AS_COMPILED_LIBRARY is ${CDT_USE_AS_COMPILED_LIBRARY}") +message(STATUS "CDT_USE_64_BIT_INDEX_TYPE is ${CDT_USE_64_BIT_INDEX_TYPE}") # Use boost for c++98 versions of c++11 containers or for Boost::rtree if(CDT_USE_BOOST) @@ -61,7 +66,10 @@ if(CDT_USE_AS_COMPILED_LIBRARY) set(cdt_headers include/CDT.h include/CDTUtils.h - include/PointRTree.h + include/KDTree.h + include/LocatorBoostRTree.h + include/LocatorNearestRandom.h + include/LocatorKDTree.h include/remove_at.hpp include/CDT.hpp include/CDTUtils.hpp @@ -114,6 +122,7 @@ target_compile_definitions( ${cdt_scope} $<$:CDT_USE_BOOST> $<$:CDT_USE_AS_COMPILED_LIBRARY> + $<$:CDT_USE_64_BIT_INDEX_TYPE> ) if(CDT_USE_BOOST) diff --git a/CDT/extras/VerifyTopology.h b/CDT/extras/VerifyTopology.h index a1ea6363..08f1f6fe 100644 --- a/CDT/extras/VerifyTopology.h +++ b/CDT/extras/VerifyTopology.h @@ -27,15 +27,15 @@ namespace CDT * * @tparam T type of vertex coordinates (e.g., float, double) */ -template -inline bool verifyTopology(const CDT::Triangulation& cdt) +template +inline bool verifyTopology(const CDT::Triangulation& cdt) { // Check if vertices' adjacent triangles contain vertex for(VertInd iV(0); iV < VertInd(cdt.vertices.size()); ++iV) { - const Vertex& v = cdt.vertices[iV]; + const TriIndVec& vTris = cdt.vertTris[iV]; typedef TriIndVec::const_iterator TriIndCit; - for(TriIndCit it = v.triangles.begin(); it != v.triangles.end(); ++it) + for(TriIndCit it = vTris.begin(); it != vTris.end(); ++it) { const array& vv = cdt.triangles[*it].vertices; if(std::find(vv.begin(), vv.end(), iV) == vv.end()) @@ -63,7 +63,7 @@ inline bool verifyTopology(const CDT::Triangulation& cdt) typedef VerticesArr3::const_iterator VCit; for(VCit it = t.vertices.begin(); it != t.vertices.end(); ++it) { - const std::vector& tt = cdt.vertices[*it].triangles; + const TriIndVec& tt = cdt.vertTris[*it]; if(std::find(tt.begin(), tt.end(), iT) == tt.end()) return false; } diff --git a/CDT/include/CDT.h b/CDT/include/CDT.h index 42ddfa5d..6dc48820 100644 --- a/CDT/include/CDT.h +++ b/CDT/include/CDT.h @@ -11,12 +11,9 @@ #define CDT_lNrmUayWQaIR5fxnsg9B #include "CDTUtils.h" -#ifdef CDT_USE_BOOST -#include "PointRTree.h" -#endif - #include "remove_at.hpp" +#include #include #include #include @@ -29,8 +26,13 @@ namespace CDT { -/// Enum of strategies for finding closest point to the newly inserted one -struct CDT_EXPORT FindingClosestPoint +/** + * Enum of strategies specifying order in which a range of vertices is inserted + * @note @ref 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 @@ -38,10 +40,8 @@ struct CDT_EXPORT FindingClosestPoint */ enum Enum { -#ifdef CDT_USE_BOOST - BoostRTree, ///< use boost::geometry::rtree -#endif - ClosestRandom, ///< pick closest from few randomly selected candidates + Randomized, ///< vertices will be inserted in random order + AsProvided, ///< vertices will be inserted in the same order as provided }; }; @@ -60,11 +60,12 @@ struct CDT_EXPORT SuperGeometryType }; /// Constant representing no valid neighbor for a triangle -const static TriInd noNeighbor(std::numeric_limits::max()); +const static TriInd noNeighbor(std::numeric_limits::max()); /// Constant representing no valid vertex for a triangle -const static VertInd noVertex(std::numeric_limits::max()); +const static VertInd noVertex(std::numeric_limits::max()); -/** type used for storing layer depths for triangles +/** + * 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. */ @@ -75,15 +76,20 @@ typedef LayerDepth BoundaryOverlapCount; * 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 +template class CDT_EXPORT Triangulation { public: - typedef std::vector > VertexVec; ///< Vertices vector - VertexVec vertices; ///< triangulation's vertices - TriangleVec triangles; ///< triangulation's triangles - EdgeUSet fixedEdges; ///< triangulation's constraints (fixed edges) + typedef std::vector > V2dVec; ///< Vertices vector + typedef std::vector VerticesTriangles; ///< Triangles by vertex + V2dVec vertices; ///< triangulation's vertices + TriangleVec triangles; ///< triangulation's triangles + EdgeUSet fixedEdges; ///< triangulation's constraints (fixed edges) + VerticesTriangles vertTris; ///< triangles adjacent to each vertex /** Stores count of overlapping boundaries for a fixed edge. If no entry is * present for an edge: no boundaries overlap. @@ -95,10 +101,22 @@ class CDT_EXPORT Triangulation unordered_map overlapCount; /*____ API _____*/ - /// Constructor + /// 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 nearPtLocator class providing locating near point for efficiently + * inserting new points + */ Triangulation( - const FindingClosestPoint::Enum closestPtMode, - const size_t nRandSamples = 10); + VertexInsertionOrder::Enum vertexInsertionOrder, + const TNearPointLocator& nearPtLocator); /** * Insert custom point-types specified by iterator range and X/Y-getters * @tparam TVertexIter iterator that dereferences to custom point type @@ -180,7 +198,8 @@ class CDT_EXPORT Triangulation private: /*____ Detail __*/ void addSuperTriangle(const Box2d& box); - void insertVertex(const V2d& pos); + void addNexVertex(const V2d& pos, const TriIndVec& tris); + void insertVertex(const VertInd iVert); void insertEdge(Edge edge); tuple intersectedTriangle( const VertInd iA, @@ -188,19 +207,13 @@ class CDT_EXPORT Triangulation const V2d& a, const V2d& b) const; /// Returns indices of three resulting triangles - std::stack - insertPointInTriangle(const V2d& pos, const TriInd iT); + std::stack insertPointInTriangle(const VertInd v, const TriInd iT); /// Returns indices of four resulting triangles std::stack - insertPointOnEdge(const V2d& pos, const TriInd iT1, const TriInd iT2); + insertPointOnEdge(const VertInd v, const TriInd iT1, const TriInd iT2); array trianglesAt(const V2d& pos) const; array walkingSearchTrianglesAt(const V2d& pos) const; TriInd walkTriangles(const VertInd startVertex, const V2d& pos) const; - VertInd - nearestVertexRand(const V2d& pos, const std::size_t nSamples) const; -#ifdef CDT_USE_BOOST - VertInd nearestVertexRtree(const V2d& pos) const; -#endif bool isFlipNeeded( const V2d& pos, const TriInd iT, @@ -217,6 +230,17 @@ class CDT_EXPORT Triangulation const VertInd iVedge2, const TriInd newNeighbor); void addAdjacentTriangle(const VertInd iVertex, const TriInd iTriangle); + void addAdjacentTriangles( + const VertInd iVertex, + const TriInd iT1, + const TriInd iT2, + const TriInd iT3); + void addAdjacentTriangles( + const VertInd iVertex, + const TriInd iT1, + const TriInd iT2, + const TriInd iT3, + const TriInd iT4); void removeAdjacentTriangle(const VertInd iVertex, const TriInd iTriangle); TriInd triangulatePseudopolygon( const VertInd ia, @@ -238,13 +262,10 @@ class CDT_EXPORT Triangulation void fixEdge(const Edge& edge); std::vector m_dummyTris; -#ifdef CDT_USE_BOOST - PointRTree m_rtree; -#endif - std::size_t m_nRandSamples; - FindingClosestPoint::Enum m_closestPtMode; + TNearPointLocator m_nearPtLocator; std::size_t m_nTargetVerts; SuperGeometryType::Enum m_superGeomType; + VertexInsertionOrder::Enum m_vertexInsertionOrder; }; /** @@ -369,15 +390,13 @@ CDT_EXPORT DuplicatesInfo RemoveDuplicatesAndRemapEdges( * - 2 for triangles in hole * - 3 for triangles in island and so on... * - * @tparam T type of vertex coordinates (e.g., float, double) - * @param vertices vertices of triangulation + * @param seed seed triangle to begin depth-peeling from * @param triangles triangles of triangulation * @param fixedEdges constraint edges of triangulation * @return vector where element at index i stores depth of i-th triangle */ -template CDT_EXPORT std::vector CalculateTriangleDepths( - const std::vector >& vertices, + const TriInd seed, const TriangleVec& triangles, const EdgeUSet& fixedEdges); @@ -393,17 +412,15 @@ CDT_EXPORT std::vector CalculateTriangleDepths( * - 2 for triangles in hole * - 3 for triangles in island and so on... * - * @tparam T type of vertex coordinates (e.g., float, double) - * @param vertices vertices of triangulation + * @param seed seed triangle to begin depth-peeling from * @param triangles triangles of triangulation * @param fixedEdges constraint edges of triangulation * @param overlapCount counts of boundary overlaps at fixed edges (map has * entries only for edges that represent overlapping boundaries) * @return vector where element at index i stores depth of i-th triangle */ -template CDT_EXPORT std::vector CalculateTriangleDepths( - const std::vector >& vertices, + const TriInd seed, const TriangleVec& triangles, const EdgeUSet& fixedEdges, const unordered_map& overlapCount); @@ -484,33 +501,62 @@ struct hash > namespace CDT { +inline size_t randomCDT(const size_t i) +{ + static std::mt19937 g(9001); + return g() % i; +} + //----------------------- // Triangulation methods //----------------------- -template +template template < typename TVertexIter, typename TGetVertexCoordX, typename TGetVertexCoordY> -void Triangulation::insertVertices( - TVertexIter first, +void Triangulation::insertVertices( + const TVertexIter first, const TVertexIter last, TGetVertexCoordX getX, TGetVertexCoordY getY) { if(vertices.empty()) + { addSuperTriangle(envelopBox(first, last, getX, getY)); - vertices.reserve(vertices.size() + std::distance(first, last)); - for(; first != last; ++first) - insertVertex(V2d::make(getX(*first), getY(*first))); + } + + const std::size_t nExistingVerts = vertices.size(); + + vertices.reserve(nExistingVerts + std::distance(first, last)); + for(TVertexIter it = first; it != last; ++it) + addNexVertex(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 = nExistingVerts; + for(Iter it = ii.begin(); it != ii.end(); ++it, ++value) + *it = value; + std::random_shuffle(ii.begin(), ii.end(), randomCDT); + for(Iter it = ii.begin(); it != ii.end(); ++it) + insertVertex(*it); + break; + } } -template +template template < typename TEdgeIter, typename TGetEdgeVertexStart, typename TGetEdgeVertexEnd> -void Triangulation::insertEdges( +void Triangulation::insertEdges( TEdgeIter first, const TEdgeIter last, TGetEdgeVertexStart getStart, diff --git a/CDT/include/CDT.hpp b/CDT/include/CDT.hpp index fd12e091..1b372075 100644 --- a/CDT/include/CDT.hpp +++ b/CDT/include/CDT.hpp @@ -33,18 +33,33 @@ array arr3(const T& v0, const T& v1, const T& v2) } // namespace detail -template -Triangulation::Triangulation( - const FindingClosestPoint::Enum closestPtMode, - const size_t nRandSamples) - : m_nRandSamples(nRandSamples) - , m_closestPtMode(closestPtMode) - , m_nTargetVerts(0) +template +Triangulation::Triangulation() + : m_nTargetVerts(0) , m_superGeomType(SuperGeometryType::SuperTriangle) + , m_vertexInsertionOrder(VertexInsertionOrder::Randomized) {} -template -void Triangulation::changeNeighbor( +template +Triangulation::Triangulation( + const VertexInsertionOrder::Enum vertexInsertionOrder) + : m_nTargetVerts(0) + , m_superGeomType(SuperGeometryType::SuperTriangle) + , m_vertexInsertionOrder(vertexInsertionOrder) +{} + +template +Triangulation::Triangulation( + VertexInsertionOrder::Enum vertexInsertionOrder, + const TNearPointLocator& nearPtLocator) + : m_nTargetVerts(0) + , m_nearPtLocator(nearPtLocator) + , m_superGeomType(SuperGeometryType::SuperTriangle) + , m_vertexInsertionOrder(vertexInsertionOrder) +{} + +template +void Triangulation::changeNeighbor( const TriInd iT, const VertInd iVedge1, const VertInd iVedge2, @@ -54,8 +69,8 @@ void Triangulation::changeNeighbor( t.neighbors[opposedTriangleInd(t, iVedge1, iVedge2)] = newNeighbor; } -template -void Triangulation::eraseDummies() +template +void Triangulation::eraseDummies() { if(m_dummyTris.empty()) return; @@ -73,11 +88,10 @@ void Triangulation::eraseDummies() triangles.erase(triangles.end() - dummySet.size(), triangles.end()); // remap adjacent triangle indices for vertices - typedef typename VertexVec::iterator VertIt; - for(VertIt v = vertices.begin(); v != vertices.end(); ++v) + typedef typename VerticesTriangles::iterator VertTrisIt; + for(VertTrisIt vTris = vertTris.begin(); vTris != vertTris.end(); ++vTris) { - TriIndVec& vTris = v->triangles; - for(TriIndVec::iterator iT = vTris.begin(); iT != vTris.end(); ++iT) + for(TriIndVec::iterator iT = vTris->begin(); iT != vTris->end(); ++iT) *iT = triIndMap[*iT]; } // remap neighbor indices for triangles @@ -91,8 +105,8 @@ void Triangulation::eraseDummies() m_dummyTris = std::vector(); } -template -void Triangulation::eraseSuperTriangleVertices() +template +void Triangulation::eraseSuperTriangleVertices() { if(m_superGeomType != SuperGeometryType::SuperTriangle) return; @@ -109,11 +123,12 @@ void Triangulation::eraseSuperTriangleVertices() } fixedEdges = updatedFixedEdges; - vertices = std::vector >(vertices.begin() + 3, vertices.end()); + vertices = std::vector >(vertices.begin() + 3, vertices.end()); + vertTris = VerticesTriangles(vertTris.begin() + 3, vertTris.end()); } -template -void Triangulation::eraseSuperTriangle() +template +void Triangulation::eraseSuperTriangle() { if(m_superGeomType != SuperGeometryType::SuperTriangle) return; @@ -128,22 +143,21 @@ void Triangulation::eraseSuperTriangle() eraseSuperTriangleVertices(); } -template -void Triangulation::eraseOuterTriangles() +template +void Triangulation::eraseOuterTriangles() { // make dummy triangles adjacent to super-triangle's vertices - const std::stack seed( - std::deque(1, vertices[0].triangles.front())); + const std::stack seed(std::deque(1, vertTris[0].front())); const TriIndUSet toErase = growToBoundary(seed); eraseTrianglesAtIndices(toErase.begin(), toErase.end()); eraseSuperTriangleVertices(); } -template -void Triangulation::eraseOuterTrianglesAndHoles() +template +void Triangulation::eraseOuterTrianglesAndHoles() { - const std::vector triDepths = - CalculateTriangleDepths(vertices, triangles, fixedEdges, overlapCount); + const std::vector triDepths = CalculateTriangleDepths( + vertTris[0].front(), triangles, fixedEdges, overlapCount); TriIndVec toErase; toErase.reserve(triangles.size()); @@ -157,9 +171,9 @@ void Triangulation::eraseOuterTrianglesAndHoles() eraseSuperTriangleVertices(); } -template +template template -void Triangulation::eraseTrianglesAtIndices( +void Triangulation::eraseTrianglesAtIndices( TriIndexIter first, TriIndexIter last) { @@ -168,21 +182,20 @@ void Triangulation::eraseTrianglesAtIndices( eraseDummies(); } -template -void Triangulation::initializedWithCustomSuperGeometry() +template +void Triangulation::initializedWithCustomSuperGeometry() { -#ifdef CDT_USE_BOOST for(std::size_t i = 0; i < vertices.size(); ++i) { - m_rtree.addPoint(vertices[i].pos, VertInd(i)); + m_nearPtLocator.addPoint(VertInd(i), vertices); } -#endif m_nTargetVerts = vertices.size(); m_superGeomType = SuperGeometryType::Custom; } -template -TriIndUSet Triangulation::growToBoundary(std::stack seeds) const +template +TriIndUSet Triangulation::growToBoundary( + std::stack seeds) const { TriIndUSet traversed; while(!seeds.empty()) @@ -204,8 +217,8 @@ TriIndUSet Triangulation::growToBoundary(std::stack seeds) const return traversed; } -template -void Triangulation::makeDummy(const TriInd iT) +template +void Triangulation::makeDummy(const TriInd iT) { const Triangle& t = triangles[iT]; @@ -220,8 +233,8 @@ void Triangulation::makeDummy(const TriInd iT) m_dummyTris.push_back(iT); } -template -TriInd Triangulation::addTriangle(const Triangle& t) +template +TriInd Triangulation::addTriangle(const Triangle& t) { if(m_dummyTris.empty()) { @@ -234,8 +247,8 @@ TriInd Triangulation::addTriangle(const Triangle& t) return nxtDummy; } -template -TriInd Triangulation::addTriangle() +template +TriInd Triangulation::addTriangle() { if(m_dummyTris.empty()) { @@ -250,14 +263,15 @@ TriInd Triangulation::addTriangle() return nxtDummy; } -template -void Triangulation::insertEdges(const std::vector& edges) +template +void Triangulation::insertEdges( + const std::vector& edges) { insertEdges(edges.begin(), edges.end(), edge_get_v1, edge_get_v2); } -template -void Triangulation::fixEdge(const Edge& edge) +template +void Triangulation::fixEdge(const Edge& edge) { if(!fixedEdges.insert(edge).second) { @@ -265,24 +279,25 @@ void Triangulation::fixEdge(const Edge& edge) } } -template -void Triangulation::insertEdge(Edge edge) +template +void Triangulation::insertEdge(Edge edge) { const VertInd iA = edge.v1(); VertInd iB = edge.v2(); if(iA == iB) // edge connects a vertex to itself return; - const Vertex& a = vertices[iA]; - const Vertex& b = vertices[iB]; - if(verticesShareEdge(a, b)) + 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(iA, iB)); return; } TriInd iT; VertInd iVleft, iVright; - tie(iT, iVleft, iVright) = - intersectedTriangle(iA, a.triangles, a.pos, b.pos); + tie(iT, iVleft, iVright) = intersectedTriangle(iA, aTris, a, b); // if one of the triangle vertices is on the edge, move edge start if(iT == noNeighbor) { @@ -300,13 +315,13 @@ void Triangulation::insertEdge(Edge edge) const TriInd iTopo = opposedTriangle(t, iV); const Triangle& tOpo = triangles[iTopo]; const VertInd iVopo = opposedVertex(tOpo, iT); - const Vertex vOpo = vertices[iVopo]; + const V2d vOpo = vertices[iVopo]; intersected.push_back(iTopo); iT = iTopo; t = triangles[iT]; - PtLineLocation::Enum loc = locatePointLine(vOpo.pos, a.pos, b.pos); + const PtLineLocation::Enum loc = locatePointLine(vOpo, a, b); if(loc == PtLineLocation::Left) { ptsLeft.push_back(iVopo); @@ -348,8 +363,9 @@ void Triangulation::insertEdge(Edge edge) * - index of point on the line * - index of point on the right of the line */ -template -tuple Triangulation::intersectedTriangle( +template +tuple +Triangulation::intersectedTriangle( const VertInd iA, const std::vector& candidates, const V2d& a, @@ -363,10 +379,8 @@ tuple Triangulation::intersectedTriangle( const Index i = vertexInd(t, iA); const VertInd iP1 = t.vertices[cw(i)]; const VertInd iP2 = t.vertices[ccw(i)]; - const PtLineLocation::Enum locP1 = - locatePointLine(vertices[iP1].pos, a, b); - const PtLineLocation::Enum locP2 = - locatePointLine(vertices[iP2].pos, a, b); + const PtLineLocation::Enum locP1 = locatePointLine(vertices[iP1], a, b); + const PtLineLocation::Enum locP2 = locatePointLine(vertices[iP2], a, b); if(locP2 == PtLineLocation::Right) { if(locP1 == PtLineLocation::OnLine) @@ -379,8 +393,8 @@ tuple Triangulation::intersectedTriangle( "edge. Note: can be caused by duplicate points."); } -template -void Triangulation::addSuperTriangle(const Box2d& box) +template +void Triangulation::addSuperTriangle(const Box2d& box) { m_nTargetVerts = 3; m_superGeomType = SuperGeometryType::SuperTriangle; @@ -396,31 +410,37 @@ void Triangulation::addSuperTriangle(const Box2d& box) 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}; - vertices.push_back(Vertex::make(posV1, TriInd(0))); - vertices.push_back(Vertex::make(posV2, TriInd(0))); - vertices.push_back(Vertex::make(posV3, TriInd(0))); + addNexVertex(posV1, TriIndVec(1, TriInd(0))); + addNexVertex(posV2, TriIndVec(1, TriInd(0))); + addNexVertex(posV3, TriIndVec(1, TriInd(0))); const Triangle superTri = { {VertInd(0), VertInd(1), VertInd(2)}, {noNeighbor, noNeighbor, noNeighbor}}; addTriangle(superTri); -#ifdef CDT_USE_BOOST - if(m_closestPtMode == FindingClosestPoint::BoostRTree) - { - m_rtree.addPoint(posV1, VertInd(0)); - m_rtree.addPoint(posV2, VertInd(1)); - m_rtree.addPoint(posV3, VertInd(2)); - } -#endif + + m_nearPtLocator.addPoint(VertInd(0), vertices); + m_nearPtLocator.addPoint(VertInd(1), vertices); + m_nearPtLocator.addPoint(VertInd(2), vertices); } -template -void Triangulation::insertVertex(const V2d& pos) +template +void Triangulation::addNexVertex( + const V2d& pos, + const TriIndVec& tris) { - const VertInd iVert(vertices.size()); - array trisAt = walkingSearchTrianglesAt(pos); + vertices.push_back(pos); + vertTris.push_back(tris); +} + +template +void Triangulation::insertVertex(const VertInd iVert) +{ + const V2d& v = vertices[iVert]; + array trisAt = walkingSearchTrianglesAt(v); std::stack triStack = - trisAt[1] == noNeighbor ? insertPointInTriangle(pos, trisAt[0]) - : insertPointOnEdge(pos, trisAt[0], trisAt[1]); + trisAt[1] == noNeighbor + ? insertPointInTriangle(iVert, trisAt[0]) + : insertPointOnEdge(iVert, trisAt[0], trisAt[1]); while(!triStack.empty()) { const TriInd iT = triStack.top(); @@ -430,17 +450,15 @@ void Triangulation::insertVertex(const V2d& pos) const TriInd iTopo = opposedTriangle(t, iVert); if(iTopo == noNeighbor) continue; - if(isFlipNeeded(pos, iT, iTopo, iVert)) + if(isFlipNeeded(v, iT, iTopo, iVert)) { flipEdge(iT, iTopo); triStack.push(iT); triStack.push(iTopo); } } -#ifdef CDT_USE_BOOST - if(m_closestPtMode == FindingClosestPoint::BoostRTree) - m_rtree.addPoint(pos, iVert); -#endif + + m_nearPtLocator.addPoint(iVert, vertices); } /*! @@ -453,8 +471,8 @@ void Triangulation::insertVertex(const V2d& pos) * as the non-super-tri shared vertex * 3. None of the vertices are super-tri: normal circumcircle test */ -template -bool Triangulation::isFlipNeeded( +template +bool Triangulation::isFlipNeeded( const V2d& pos, const TriInd iT, const TriInd iTopo, @@ -465,13 +483,9 @@ bool Triangulation::isFlipNeeded( const VertInd iVopo = tOpo.vertices[i]; const VertInd iVcw = tOpo.vertices[cw(i)]; const VertInd iVccw = tOpo.vertices[ccw(i)]; - - if(fixedEdges.count(Edge(iVcw, iVccw))) - return false; // flip not needed if the original edge is fixed - - const V2d& v1 = vertices[iVcw].pos; - const V2d& v2 = vertices[iVopo].pos; - const V2d& v3 = vertices[iVccw].pos; + const V2d& v1 = vertices[iVcw]; + const V2d& v2 = vertices[iVopo]; + const V2d& v3 = vertices[iVccw]; if(m_superGeomType == SuperGeometryType::SuperTriangle) { if(iVert < 3 || iVopo < 3) // flip-candidate edge touches super-triangle @@ -508,11 +522,11 @@ bool Triangulation::isFlipNeeded( * v1 ___________________ v2 * n1 */ -template -std::stack -Triangulation::insertPointInTriangle(const V2d& pos, const TriInd iT) +template +std::stack Triangulation::insertPointInTriangle( + const VertInd v, + const TriInd iT) { - const VertInd v(vertices.size()); const TriInd iNewT1 = addTriangle(); const TriInd iNewT2 = addTriangle(); @@ -527,7 +541,7 @@ Triangulation::insertPointInTriangle(const V2d& pos, const TriInd 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 - vertices.push_back(Vertex::makeInTriangle(pos, iT, iNewT1, iNewT2)); + addAdjacentTriangles(v, iT, iNewT1, iNewT2); // adjust lists of adjacent triangles for v1, v2, v3 addAdjacentTriangle(v1, iNewT2); addAdjacentTriangle(v2, iNewT1); @@ -558,13 +572,12 @@ Triangulation::insertPointInTriangle(const V2d& pos, const TriInd iT) * \|/ * T2 (bottom) v3 */ -template -std::stack Triangulation::insertPointOnEdge( - const V2d& pos, +template +std::stack Triangulation::insertPointOnEdge( + const VertInd v, const TriInd iT1, const TriInd iT2) { - const VertInd v(vertices.size()); const TriInd iTnew1 = addTriangle(); const TriInd iTnew2 = addTriangle(); @@ -587,7 +600,7 @@ std::stack Triangulation::insertPointOnEdge( 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 - vertices.push_back(Vertex::makeOnEdge(pos, iT1, iTnew2, iT2, iTnew1)); + addAdjacentTriangles(v, iT1, iTnew2, iT2, iTnew1); // adjust neighboring triangles and vertices changeNeighbor(n4, iT1, iTnew1); changeNeighbor(n2, iT2, iTnew2); @@ -606,16 +619,17 @@ std::stack Triangulation::insertPointOnEdge( return newTriangles; } -template -array Triangulation::trianglesAt(const V2d& pos) const +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]].pos; - const V2d v2 = vertices[t.vertices[1]].pos; - const V2d v3 = vertices[t.vertices[2]].pos; + 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; @@ -627,13 +641,13 @@ array Triangulation::trianglesAt(const V2d& pos) const throw std::runtime_error("No triangle was found at position"); } -template -TriInd Triangulation::walkTriangles( +template +TriInd Triangulation::walkTriangles( const VertInd startVertex, const V2d& pos) const { // begin walk in search of triangle at pos - TriInd currTri = vertices[startVertex].triangles[0]; + TriInd currTri = vertTris[startVertex][0]; #ifdef CDT_USE_BOOST TriIndFlatUSet visited; #else @@ -649,8 +663,8 @@ TriInd Triangulation::walkTriangles( for(Index i_(0); i_ < Index(3); ++i_) { const Index i((i_ + offset) % 3); - const V2d vStart = vertices[t.vertices[i]].pos; - const V2d vEnd = vertices[t.vertices[ccw(i)]].pos; + 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 && @@ -666,26 +680,19 @@ TriInd Triangulation::walkTriangles( return currTri; } -template -array -Triangulation::walkingSearchTrianglesAt(const V2d& pos) const +template +array Triangulation::walkingSearchTrianglesAt( + const V2d& pos) const { array out = {noNeighbor, noNeighbor}; - // Query RTree for a vertex close to pos, to start the search -#ifdef CDT_USE_BOOST - const VertInd startVertex = - m_closestPtMode == FindingClosestPoint::BoostRTree - ? nearestVertexRtree(pos) - : nearestVertexRand(pos, m_nRandSamples); -#else - const VertInd startVertex = nearestVertexRand(pos, m_nRandSamples); -#endif + // 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]].pos; - const V2d v2 = vertices[t.vertices[1]].pos; - const V2d v3 = vertices[t.vertices[2]].pos; + 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"); @@ -695,35 +702,6 @@ Triangulation::walkingSearchTrianglesAt(const V2d& pos) const return out; } -#ifdef CDT_USE_BOOST -template -VertInd Triangulation::nearestVertexRtree(const V2d& pos) const -{ - return m_rtree.nearestPoint(pos); -} -#endif - -template -VertInd Triangulation::nearestVertexRand( - const V2d& pos, - const std::size_t nSamples) const -{ - // start search at a vertex close to pos based on random sampling - VertInd out(detail::randGen() % vertices.size()); - T minDist = distance(vertices[out].pos, pos); - for(std::size_t iSample = 0; iSample < nSamples; ++iSample) - { - const VertInd candidate(detail::randGen() % vertices.size()); - const T candidateDist = distance(vertices[candidate].pos, pos); - if(candidateDist < minDist) - { - minDist = candidateDist; - out = candidate; - } - } - return out; -} - /* Flip edge between T and Topo: * * v4 | - old edge @@ -740,8 +718,10 @@ VertInd Triangulation::nearestVertexRand( * \|/ * v2 */ -template -void Triangulation::flipEdge(const TriInd iT, const TriInd iTopo) +template +void Triangulation::flipEdge( + const TriInd iT, + const TriInd iTopo) { Triangle& t = triangles[iT]; Triangle& tOpo = triangles[iTopo]; @@ -773,8 +753,8 @@ void Triangulation::flipEdge(const TriInd iT, const TriInd iTopo) removeAdjacentTriangle(v4, iTopo); } -template -void Triangulation::changeNeighbor( +template +void Triangulation::changeNeighbor( const TriInd iT, const TriInd oldNeighbor, const TriInd newNeighbor) @@ -785,20 +765,50 @@ void Triangulation::changeNeighbor( t.neighbors[neighborInd(t, oldNeighbor)] = newNeighbor; } -template -void Triangulation::addAdjacentTriangle( +template +void Triangulation::addAdjacentTriangle( const VertInd iVertex, const TriInd iTriangle) { - vertices[iVertex].triangles.push_back(iTriangle); + vertTris[iVertex].push_back(iTriangle); } -template -void Triangulation::removeAdjacentTriangle( +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 = vertices[iVertex].triangles; + std::vector& tris = vertTris[iVertex]; tris.erase(std::find(tris.begin(), tris.end(), iTriangle)); } @@ -815,8 +825,8 @@ splitPseudopolygon(const VertInd vi, const std::vector& points) return out; } -template -TriInd Triangulation::triangulatePseudopolygon( +template +TriInd Triangulation::triangulatePseudopolygon( const VertInd ia, const VertInd ib, const std::vector& points) @@ -854,36 +864,36 @@ TriInd Triangulation::triangulatePseudopolygon( return iT; } -template -VertInd Triangulation::findDelaunayPoint( +template +VertInd Triangulation::findDelaunayPoint( const VertInd ia, const VertInd ib, const std::vector& points) const { assert(!points.empty()); - const V2d& a = vertices[ia].pos; - const V2d& b = vertices[ib].pos; + const V2d& a = vertices[ia]; + const V2d& b = vertices[ib]; VertInd ic = points.front(); - V2d c = vertices[ic].pos; + V2d c = vertices[ic]; typedef std::vector::const_iterator CIt; for(CIt it = points.begin() + 1; it != points.end(); ++it) { - const V2d v = vertices[*it].pos; + const V2d v = vertices[*it]; if(!isInCircumcircle(v, a, b, c)) continue; ic = *it; - c = vertices[ic].pos; + c = vertices[ic]; } return ic; } -template -TriInd Triangulation::pseudopolyOuterTriangle( +template +TriInd Triangulation::pseudopolyOuterTriangle( const VertInd ia, const VertInd ib) const { - const std::vector& aTris = vertices[ia].triangles; - const std::vector& bTris = vertices[ib].triangles; + 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()) @@ -891,8 +901,9 @@ TriInd Triangulation::pseudopolyOuterTriangle( return noNeighbor; } -template -void Triangulation::insertVertices(const std::vector >& newVertices) +template +void Triangulation::insertVertices( + const std::vector >& newVertices) { return insertVertices( newVertices.begin(), newVertices.end(), getX_V2d, getY_V2d); @@ -997,16 +1008,15 @@ TriIndUSet PeelLayer( return behindBoundary; } -template std::vector CalculateTriangleDepths( - const std::vector >& vertices, + const TriInd seed, const TriangleVec& triangles, const EdgeUSet& fixedEdges, const unordered_map& overlapCount) { std::vector triDepths( triangles.size(), std::numeric_limits::max()); - std::stack seeds(TriDeque(1, vertices[0].triangles.front())); + std::stack seeds(TriDeque(1, seed)); LayerDepth layerDepth = 0; LayerDepth deepestSeedDepth = 0; @@ -1032,15 +1042,14 @@ std::vector CalculateTriangleDepths( return triDepths; } -template std::vector CalculateTriangleDepths( - const std::vector >& vertices, + const TriInd seed, const TriangleVec& triangles, const EdgeUSet& fixedEdges) { std::vector triDepths( triangles.size(), std::numeric_limits::max()); - std::stack seeds(TriDeque(1, vertices[0].triangles.front())); + std::stack seeds(TriDeque(1, seed)); LayerDepth layerDepth = 0; do diff --git a/CDT/include/CDTUtils.h b/CDT/include/CDTUtils.h index 3a153be9..431d4cf9 100644 --- a/CDT/include/CDTUtils.h +++ b/CDT/include/CDTUtils.h @@ -123,20 +123,26 @@ CDT_EXPORT 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 std::size_t VertInd; +typedef IndexSizeType VertInd; /// Triangle index -typedef std::size_t TriInd; +typedef IndexSizeType TriInd; #else /// Index in triangle BOOST_STRONG_TYPEDEF(unsigned char, Index); /// Vertex index -BOOST_STRONG_TYPEDEF(std::size_t, VertInd); +BOOST_STRONG_TYPEDEF(IndexSizeType, VertInd); /// Triangle index -BOOST_STRONG_TYPEDEF(std::size_t, TriInd); +BOOST_STRONG_TYPEDEF(IndexSizeType, TriInd); #endif typedef std::vector TriIndVec; ///< Vector of triangle indices @@ -154,7 +160,7 @@ struct CDT_EXPORT Box2d /// Bounding box of a collection of custom 2D points given coordinate getters template < typename T, - typename TVertexIter, + typename TVertexIter, typename TGetVertexCoordX, typename TGetVertexCoordY> CDT_EXPORT Box2d envelopBox( @@ -179,30 +185,6 @@ CDT_EXPORT Box2d envelopBox( template CDT_EXPORT Box2d envelopBox(const std::vector >& vertices); -/// Triangulation vertex -template -struct CDT_EXPORT Vertex -{ - V2d pos; ///< vertex position - TriIndVec triangles; ///< adjacent triangles - - /// Create vertex - static Vertex make(const V2d& pos, const TriInd iTriangle); - /// Create vertex in a triangle - static Vertex makeInTriangle( - const V2d& pos, - const TriInd iT1, - const TriInd iT2, - const TriInd iT3); - /// Create vertex on an edge - static Vertex makeOnEdge( - const V2d& pos, - const TriInd iT1, - const TriInd iT2, - const TriInd iT3, - const TriInd iT4); -}; - /// Edge connecting two vertices: vertex with smaller index is always first /// \note: hash Edge is specialized at the bottom struct CDT_EXPORT Edge @@ -366,13 +348,17 @@ CDT_EXPORT bool isInCircumcircle( const V2d& v3); /// Test if two vertices share at least one common triangle -template -CDT_EXPORT bool verticesShareEdge(const Vertex& a, const Vertex& b); +CDT_EXPORT inline 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 diff --git a/CDT/include/CDTUtils.hpp b/CDT/include/CDTUtils.hpp index 7efd7220..580e41f0 100644 --- a/CDT/include/CDTUtils.hpp +++ b/CDT/include/CDTUtils.hpp @@ -38,52 +38,6 @@ Box2d envelopBox(const std::vector >& vertices) vertices.begin(), vertices.end(), getX_V2d, getY_V2d); } -//***************************************************************************** -// Vertex -//***************************************************************************** -template -Vertex Vertex::make(const V2d& pos, const TriInd iTriangle) -{ - Vertex out = {pos, std::vector(1, iTriangle)}; - return out; -} - -template -Vertex Vertex::makeInTriangle( - const V2d& pos, - const TriInd iT1, - const TriInd iT2, - const TriInd iT3) -{ - Vertex out; - out.pos = pos; - TriIndVec& vTris = out.triangles; - vTris.reserve(3); - vTris.push_back(iT1); - vTris.push_back(iT2); - vTris.push_back(iT3); - return out; -} - -template -Vertex Vertex::makeOnEdge( - const V2d& pos, - const TriInd iT1, - const TriInd iT2, - const TriInd iT3, - const TriInd iT4) -{ - Vertex out; - out.pos = pos; - TriIndVec& vTris = out.triangles; - vTris.reserve(4); - vTris.push_back(iT1); - vTris.push_back(iT2); - vTris.push_back(iT3); - vTris.push_back(iT4); - return out; -} - //***************************************************************************** // Edge //***************************************************************************** @@ -272,23 +226,39 @@ bool isInCircumcircle( return incircle(v1.x, v1.y, v2.x, v2.y, v3.x, v3.y, p.x, p.y) > T(0); } -template -bool verticesShareEdge(const Vertex& a, const Vertex& b) +CDT_INLINE_IF_HEADER_ONLY +bool verticesShareEdge(const TriIndVec& aTris, const TriIndVec& bTris) { - const std::vector& aTris = a.triangles; - const std::vector& bTris = b.triangles; 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) { - const T dx = b.x - a.x; - const T dy = b.y - a.y; - return std::sqrt(dx * dx + dy * dy); + 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/CDT/include/KDTree.h b/CDT/include/KDTree.h new file mode 100644 index 00000000..fa9a056d --- /dev/null +++ b/CDT/include/KDTree.h @@ -0,0 +1,339 @@ +/// 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/. +/// author: Andre Fecteau + +#ifndef KDTREE_KDTREE_H +#define KDTREE_KDTREE_H + +#include "CDTUtils.h" + +#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 MaxStackDepth size of statically allocated stack for nearest query. +template < + typename TCoordType, + size_t NumVerticesInLeaf = 20, + size_t MaxStackDepth = 64> +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() + : children({0, 0}) + { + data.reserve(NumVerticesInLeaf); + } + /// 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::lowest(), + std::numeric_limits::lowest())) + , m_max(point_type::make( + std::numeric_limits::max(), + std::numeric_limits::max())) + , m_isRootBoxInitialized(false) + { + 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; + + coord_type mid; + NodeSplitDirection::Enum newDir; + 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.children = {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] = {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; + 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); + // node containing point should end up on top of the stack + if(iChild == 0) + { + m_tasksStack[++iTask] = { + n.children[1], newMin, t.max, newDir, toMidSq}; + m_tasksStack[++iTask] = { + n.children[0], t.min, newMax, newDir, toMidSq}; + } + else + { + m_tasksStack[++iTask] = { + n.children[0], t.min, newMax, newDir, toMidSq}; + m_tasksStack[++iTask] = { + 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 = 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].children = {newLeaf, m_root} + : m_nodes[newRoot].children = {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].children = {newLeaf, m_root} + : m_nodes[newRoot].children = {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 = {std::min(m_min.x, p.x), std::min(m_min.y, p.y)}; + m_max = {std::max(m_max.x, p.x), std::max(m_max.y, p.y)}; + } + 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; + }; + // allocated in class (not in the 'nearest' method) for better performance + mutable CDT::array m_tasksStack; +}; + +} // namespace KDTree + +#endif // header guard \ No newline at end of file diff --git a/CDT/include/LocatorBoostRTree.h b/CDT/include/LocatorBoostRTree.h new file mode 100644 index 00000000..7817daab --- /dev/null +++ b/CDT/include/LocatorBoostRTree.h @@ -0,0 +1,55 @@ +/* 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 boost::geometry::rtree and CDT + */ + +#ifndef CDT_o2ROLeC4NKchgja4Awtg +#define CDT_o2ROLeC4NKchgja4Awtg + +#include "CDTUtils.h" + +#include +#include +#include +#include + +namespace CDT +{ + +/// R-tree holding points +template +class LocatorBoostRTree +{ +public: + /// Add point to R-tree + void addPoint(const VertInd iPoint, const std::vector >& points) + { + const V2d& pos = points[iPoint]; + m_rtree.insert(std::make_pair(Point(pos.x, pos.y), iPoint)); + } + /// Find nearest point using R-tree + VertInd nearPoint(const V2d& pos, const std::vector >&) const + { + std::vector query; + namespace bgi = boost::geometry::index; + m_rtree.query( + bgi::nearest(Point(pos.x, pos.y), 1), std::back_inserter(query)); + assert(query.size()); + return query.front().second; + } + +private: + typedef boost::geometry::model::point + Point; + typedef std::pair Value; + boost::geometry::index::rtree > + m_rtree; +}; + +} // namespace CDT + +#endif diff --git a/CDT/include/LocatorKDTree.h b/CDT/include/LocatorKDTree.h new file mode 100644 index 00000000..da026ea6 --- /dev/null +++ b/CDT/include/LocatorKDTree.h @@ -0,0 +1,42 @@ +/* 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 +class LocatorKDTree +{ +public: + /// Add point to R-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: + KDTree::KDTree m_kdTree; +}; + +} // namespace CDT + +#endif diff --git a/CDT/include/LocatorNearestRandom.h b/CDT/include/LocatorNearestRandom.h new file mode 100644 index 00000000..01e73164 --- /dev/null +++ b/CDT/include/LocatorNearestRandom.h @@ -0,0 +1,69 @@ +/* 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 + * Locator that returns nearest random point + */ + +#ifndef CDT_qG9bX9lZ2nG4uW8xF5kC +#define CDT_qG9bX9lZ2nG4uW8xF5kC + +#include "CDTUtils.h" + +#include +#include +#include + +namespace CDT +{ + +/// Locator that returns nearest random point +template +class LocatorNearestRandom +{ +public: + LocatorNearestRandom() + : m_nSamples(10) + , m_randGen(m_randGenSeed) + {} + LocatorNearestRandom(VertInd nSamples) + : m_nSamples(nSamples) + , m_randGen(m_randGenSeed) + {} + /// Add point to R-tree + void addPoint(const VertInd iPoint, const std::vector >&) + { + m_vertices.push_back(iPoint); + } + /// Find nearest random point + VertInd + nearPoint(const V2d& pos, const std::vector >& points) const + { + VertInd nearest = 0; + T minDistSq = std::numeric_limits::max(); + for(std::size_t iSample = 0; iSample < m_nSamples; ++iSample) + { + const VertInd iRand(m_randGen() % m_vertices.size()); + const VertInd iV = m_vertices[iRand]; + const T distSq = distanceSquared(points[iV], pos); + if(distSq < minDistSq) + { + minDistSq = distSq; + nearest = iV; + } + } + return nearest; + } + +private: + std::vector m_vertices; + VertInd m_nSamples; + static const unsigned int m_randGenSeed = 9001; + mutable mt19937 m_randGen; +}; + +} // namespace CDT + +#endif diff --git a/CDT/include/PointRTree.h b/CDT/include/PointRTree.h deleted file mode 100644 index 4f86be50..00000000 --- a/CDT/include/PointRTree.h +++ /dev/null @@ -1,122 +0,0 @@ -/* 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 boost::geometry::rtree and CDT - */ - -#ifndef CDT_o2ROLeC4NKchgja4Awtg -#define CDT_o2ROLeC4NKchgja4Awtg - -#include "CDTUtils.h" - -#include -#include - -// Specialize boost::geometry point traits for CDT's 2D vector type -namespace boost -{ -namespace geometry -{ -namespace traits -{ - -// point traits -/// Point tag -template -struct tag > -{ - typedef point_tag type; ///< point tag -}; - -/// Dimension (2d) -template -struct dimension > : boost::mpl::int_<2> -{}; - -/// Coordinate type -template -struct coordinate_type > -{ - typedef T type; ///< type -}; - -/// Cartesian coordinate system -template -struct coordinate_system > -{ - typedef cs::cartesian type; ///< cartesian -}; - -// point access -/// X -template -struct access, 0> -{ - /// X-getter - static inline T get(CDT::V2d const& p) - { - return p.x; - } - /// X-setter - static inline void set(CDT::V2d& p, T const& value) - { - p.x = value; - } -}; -/// Y -template -struct access, 1> -{ - /// Y-getter - static inline T get(CDT::V2d const& p) - { - return p.y; - } - /// Y-setter - static inline void set(CDT::V2d& p, T const& value) - { - p.y = value; - } -}; - -} // namespace traits -} // namespace geometry -} // namespace boost - -namespace CDT -{ - -/// R-tree holding points -template -class PointRTree -{ -public: - /// Add point to R-tree - void addPoint(const V2d& pos, const VertInd iV) - { - m_rtree.insert(std::make_pair(pos, iV)); - } - /// Find nearest point using R-tree - VertInd nearestPoint(const V2d& pos) const - { - std::vector query; - namespace bgi = boost::geometry::index; - m_rtree.query(bgi::nearest(pos, 1), std::back_inserter(query)); - assert(query.size()); - return query.front().second; - } - -private: - typedef std::pair, VertInd> VertexPos; - boost::geometry::index:: - rtree > - m_rtree; -}; - -} // namespace CDT - -#endif diff --git a/CDT/src/CDT.cpp b/CDT/src/CDT.cpp index 1952a7f9..c461d435 100644 --- a/CDT/src/CDT.cpp +++ b/CDT/src/CDT.cpp @@ -15,18 +15,27 @@ #include "CDT.hpp" #include "CDTUtils.hpp" +#ifdef CDT_USE_BOOST +#include "LocatorBoostRTree.h" +#endif +#include "LocatorKDTree.h" +#include "LocatorNearestRandom.h" namespace CDT { -template class Triangulation; -template class Triangulation; +#ifdef CDT_USE_BOOST +template class Triangulation >; +template class Triangulation >; +#endif +template class Triangulation >; +template class Triangulation >; +template class Triangulation >; +template class Triangulation >; template struct V2d; template struct V2d; template struct Box2d; template struct Box2d; -template struct Vertex; -template struct Vertex; template Box2d envelopBox(const std::vector >&); template Box2d envelopBox(const std::vector >&); @@ -41,15 +50,6 @@ template DuplicatesInfo RemoveDuplicatesAndRemapEdges( std::vector >&, std::vector&); -template std::vector CalculateTriangleDepths( - const std::vector >& vertices, - const TriangleVec& triangles, - const EdgeUSet& fixedEdges); -template std::vector CalculateTriangleDepths( - const std::vector >& vertices, - const TriangleVec& triangles, - const EdgeUSet& fixedEdges); - } // namespace CDT #endif diff --git a/visualizer/main.cpp b/visualizer/main.cpp index 2d4fbddb..71016d3e 100644 --- a/visualizer/main.cpp +++ b/visualizer/main.cpp @@ -2,6 +2,11 @@ * 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/. */ #include "CDT.h" +// #ifdef CDT_USE_BOOST +// #include "LocatorBoostRTree.h" +// #endif +// #include "LocatorNearestRandom.h" +#include "LocatorKDTree.h" #include "VerifyTopology.h" #include @@ -21,10 +26,15 @@ #include #include -typedef float CoordType; -typedef CDT::Triangulation Triangulation; +typedef double CoordType; +// #ifdef CDT_USE_BOOST +// typedef CDT::LocatorBoostRTree NearPtLocator; +// #else +// typedef CDT::LocatorNearestRandom NearPtLocator; +// #endif +typedef CDT::LocatorKDTree NearPtLocator; +typedef CDT::Triangulation Triangulation; typedef CDT::V2d V2d; -typedef CDT::Vertex Vertex; typedef CDT::Triangle Triangle; typedef CDT::Box2d Box2d; typedef CDT::Edge Edge; @@ -36,11 +46,6 @@ class CDTWidget : public QWidget public: explicit CDTWidget(QWidget* parent = NULL) : QWidget(parent) -#ifdef CDT_USE_BOOST - , m_cdt(CDT::FindingClosestPoint::BoostRTree) -#else - , m_cdt(CDT::FindingClosestPoint::ClosestRandom, 10) -#endif , m_ptLimit(9999999) , m_edgeLimit(9999999) , m_isHidePoints(false) @@ -129,14 +134,14 @@ public slots: const CoordType stZ = -std::fmax(box.max.x - box.min.x, box.max.y - box.min.y); std::size_t counter = 0; - typedef Triangulation::VertexVec::const_iterator VCit; + typedef Triangulation::V2dVec::const_iterator VCit; for(VCit v = m_cdt.vertices.begin(); v != m_cdt.vertices.end(); ++v) { const CoordType z = !m_isRemoveOuterAndHoles && !m_isRemoveOuter && !m_isHideSuperTri && counter < 3 ? stZ : 0.0; - fout << v->pos.x << ' ' << v->pos.y << ' ' << z << "\n"; + fout << v->x << ' ' << v->y << ' ' << z << "\n"; counter++; } // Write faces @@ -180,11 +185,7 @@ public slots: void updateCDT() { -#ifdef CDT_USE_BOOST - m_cdt = Triangulation(CDT::FindingClosestPoint::BoostRTree); -#else - m_cdt = Triangulation(CDT::FindingClosestPoint::ClosestRandom, 10); -#endif + m_cdt = Triangulation(); if(!m_points.empty()) { std::vector pts = @@ -287,9 +288,9 @@ public slots: if(t->vertices[0] > 2 && t->vertices[1] > 2 && t->vertices[2] > 2) continue; - const V2d& v1 = m_cdt.vertices[t->vertices[0]].pos; - const V2d& v2 = m_cdt.vertices[t->vertices[1]].pos; - const V2d& v3 = m_cdt.vertices[t->vertices[2]].pos; + const V2d& v1 = m_cdt.vertices[t->vertices[0]]; + const V2d& v2 = m_cdt.vertices[t->vertices[1]]; + const V2d& v3 = m_cdt.vertices[t->vertices[2]]; const QPointF pt1 = sceneToScreen(v1); const QPointF pt2 = sceneToScreen(v2); const QPointF pt3 = sceneToScreen(v3); @@ -314,9 +315,9 @@ public slots: continue; } } - const V2d& v1 = m_cdt.vertices[t->vertices[0]].pos; - const V2d& v2 = m_cdt.vertices[t->vertices[1]].pos; - const V2d& v3 = m_cdt.vertices[t->vertices[2]].pos; + const V2d& v1 = m_cdt.vertices[t->vertices[0]]; + const V2d& v2 = m_cdt.vertices[t->vertices[1]]; + const V2d& v3 = m_cdt.vertices[t->vertices[2]]; const CDT::array pts = { sceneToScreen(v1), sceneToScreen(v2), sceneToScreen(v3)}; p.drawPolygon(pts.data(), pts.size()); @@ -327,8 +328,8 @@ public slots: typedef CDT::EdgeUSet::const_iterator ECit; for(ECit e = m_cdt.fixedEdges.begin(); e != m_cdt.fixedEdges.end(); ++e) { - const V2d& v1 = m_cdt.vertices[e->v1()].pos; - const V2d& v2 = m_cdt.vertices[e->v2()].pos; + const V2d& v1 = m_cdt.vertices[e->v1()]; + const V2d& v2 = m_cdt.vertices[e->v2()]; p.drawLine(sceneToScreen(v1), sceneToScreen(v2)); } @@ -340,7 +341,7 @@ public slots: p.setPen(pen); for(std::size_t i = 0; i < m_cdt.vertices.size(); ++i) { - p.drawPoint(sceneToScreen(m_cdt.vertices[i].pos)); + p.drawPoint(sceneToScreen(m_cdt.vertices[i])); } // last added point if(m_ptLimit <= m_points.size()) @@ -348,8 +349,8 @@ public slots: pen.setColor(QColor(200, 50, 50)); pen.setWidthF(9.0); p.setPen(pen); - const Vertex& v = m_cdt.vertices.back(); - p.drawPoint(sceneToScreen(v.pos)); + const V2d& v = m_cdt.vertices.back(); + p.drawPoint(sceneToScreen(v)); } }