diff --git a/CDT/include/CDTUtils.h b/CDT/include/CDTUtils.h index 89dae343..8d594e4b 100644 --- a/CDT/include/CDTUtils.h +++ b/CDT/include/CDTUtils.h @@ -44,6 +44,7 @@ typedef char couldnt_parse_cxx_standard[-1]; ///< Error: couldn't parse standard #include #include #include +#include #include #ifdef CDT_USE_STRONG_TYPING @@ -243,6 +244,7 @@ inline Edge edge_make(VertInd iV1, VertInd iV2) } typedef std::vector EdgeVec; ///< Vector of edges +typedef std::queue EdgeQue; ///< Queue 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 diff --git a/CDT/include/Triangulation.h b/CDT/include/Triangulation.h index 7f877138..dfd2a89e 100644 --- a/CDT/include/Triangulation.h +++ b/CDT/include/Triangulation.h @@ -83,6 +83,22 @@ struct CDT_EXPORT IntersectingConstraintEdges }; }; +/** + * Enum of strategies for triangles refinement + */ +struct CDT_EXPORT RefineTriangles +{ + /** + * The Enum itself + * @note needed to pre c++11 compilers that don't support 'class enum' + */ + enum Enum + { + ByAngle, ///< constraint minimum triangles angle + ByArea, ///< constraint maximum triangles area + }; +}; + /** * Type used for storing layer depths for triangles * @note LayerDepth should support 60K+ layers, which could be to much or @@ -110,8 +126,14 @@ class CDT_EXPORT Triangulation { public: typedef std::vector > V2dVec; ///< Vertices vector + typedef std::vector boolVec; ///< Steiner Vertices flag V2dVec vertices; ///< triangulation's vertices - TriangleVec triangles; ///< triangulation's triangles + boolVec isSteinerVertex; ///< triangulation's vertices Steiner point flag + IndexSizeType + maxSteinerPoints; ///< triangulation's maximum number of Steiner + IndexSizeType numOfSteinerPoints; ///< triangulation's maximum number of + ///< Steiner points to be added + TriangleVec triangles; ///< triangulation's triangles EdgeUSet fixedEdges; ///< triangulation's constraints (fixed edges) /** Stores count of overlapping boundaries for a fixed edge. If no entry is @@ -275,6 +297,16 @@ class CDT_EXPORT Triangulation * @tparam edges edges to conform to */ void conformToEdges(const std::vector& edges); + /** + * Traingles refinement by removing bad triangles + * @note bad triangles don't fulfill constraints defined by the user + * @param refinement_constrain refinement strategy that is used to identify + * bad triangles + * @param threshold threshold value for refinement + */ + void refineTriangles( + RefineTriangles::Enum refinementConstrain = RefineTriangles::ByAngle, + T threshold = 20 / 180.0 * M_PI); /** * Erase triangles adjacent to super triangle * @@ -360,7 +392,7 @@ class CDT_EXPORT Triangulation private: /*____ Detail __*/ void addSuperTriangle(const Box2d& box); - void addNewVertex(const V2d& pos, TriInd iT); + void addNewVertex(const V2d& pos, TriInd iT, bool isSteiner = false); void insertVertex(VertInd iVert); void insertVertex(VertInd iVert, VertInd walkStart); void ensureDelaunayByEdgeFlips( @@ -488,6 +520,32 @@ class CDT_EXPORT Triangulation VertInd iV2, VertInd iV3, VertInd iV4) const; + TriInd edgeTriangle(Edge edge) const; + /// Checks if edge e is encroached by vertex v + bool isEncroached(const V2d& v, Edge e) const; + bool isBadTriangle( + const Triangle& tri, + RefineTriangles::Enum refinement, + T threshold) const; + V2d circumcenter(const Triangle& tri) const; + /// Search in all fixed edges to find encroached edges, each fixed edge is + /// checked against its opposite vertices + /// Returns queue of encroached edges + EdgeQue detectEncroachedEdges(); + /// Search in all fixed edges to find encroached edges, each fixed edge is + /// checked against its opposite vertices and vertex v + /// Returns queue of encroached edges + EdgeQue detectEncroachedEdges(const V2d& v); + /// Recursively split encroached edges + /// returns vector of badly shaped triangles and number of splits + std::pair resolveEncroachedEdges( + EdgeQue encroachedEdges, + const V2d& v = {}, + bool validV = false, + bool fillBadTriangles = false, + RefineTriangles::Enum refinementConstrain = {}, + T badTriangleThreshold = {}); + VertInd splitEncroachedEdge(Edge e, TriInd iT, TriInd iTopo); void changeNeighbor(TriInd iT, TriInd oldNeighbor, TriInd newNeighbor); void changeNeighbor( TriInd iT, @@ -745,6 +803,92 @@ void Triangulation::conformToEdges( eraseDummies(); } +template +void Triangulation::refineTriangles( + RefineTriangles::Enum refinementConstrain, + T threshold) +{ + if(isFinalized()) + { + throw std::runtime_error("Triangulation was finalized with 'erase...' " + "method. Refinement is not possible"); + } + tryInitNearestPointLocator(); + resolveEncroachedEdges(detectEncroachedEdges()); + + std::queue badTriangles; + for(TriInd iT(0), n = triangles.size(); iT < n; ++iT) + { + const Triangle& t = triangles[iT]; + if(t.vertices[0] < 3 || t.vertices[1] < 3 || t.vertices[2] < 3) + continue; + + if(isBadTriangle(t, refinementConstrain, threshold)) + { + const V2d vert = circumcenter(t); + if(locatePointTriangle( + vert, vertices[0], vertices[1], vertices[2]) != + PtTriLocation::Outside) + { + badTriangles.push(iT); + } + } + } + + while(!badTriangles.empty()) + { + TriInd iT = badTriangles.front(); + const Triangle& t = triangles[iT]; + badTriangles.pop(); + if(!isBadTriangle(t, refinementConstrain, threshold) || + numOfSteinerPoints >= maxSteinerPoints) + { + continue; + } + const V2d vert = circumcenter(t); + if(locatePointTriangle(vert, vertices[0], vertices[1], vertices[2]) == + PtTriLocation::Outside) + { + continue; + } + TriIndVec badTris = resolveEncroachedEdges( + detectEncroachedEdges(vert), + vert, + true, + true, + refinementConstrain, + threshold) + .first; + + for(IndexSizeType i(0); i < TriInd(badTris.size()); ++i) + { + badTriangles.push(badTris[i]); + } + + if(badTris.empty() && numOfSteinerPoints < maxSteinerPoints) + { + const VertInd iVert = static_cast(vertices.size()); + addNewVertex(vert, noNeighbor, true); + insertVertex(iVert); + TriInd start = m_vertTris[iVert]; + TriInd currTri = start; + do + { + const Triangle& t = triangles[currTri]; + if(isBadTriangle(t, refinementConstrain, threshold)) + { + badTriangles.push(currTri); + } + currTri = t.next(iVert).first; + } while(currTri != start); + } + else + { + badTriangles.push(iT); + } + } +} + } // namespace CDT #ifndef CDT_USE_AS_COMPILED_LIBRARY diff --git a/CDT/include/Triangulation.hpp b/CDT/include/Triangulation.hpp index 3b626a2a..5b2075f1 100644 --- a/CDT/include/Triangulation.hpp +++ b/CDT/include/Triangulation.hpp @@ -426,7 +426,9 @@ namespace detail template T lerp(const T& a, const T& b, const T t) { - return (T(1) - t) * a + t * b; + if((a <= 0 && b >= 0) || (a >= 0 && b <= 0)) + return (T(1) - t) * a + t * b; + return a + t * (b - a); } // Precondition: ab and cd intersect normally @@ -1018,9 +1020,11 @@ void Triangulation::addSuperTriangle(const Box2d& box) template void Triangulation::addNewVertex( const V2d& pos, - const TriInd iT) + const TriInd iT, + const bool isSteiner) { vertices.push_back(pos); + isSteinerVertex.push_back(isSteiner); m_vertTris.push_back(iT); } @@ -1262,6 +1266,282 @@ bool Triangulation::isFlipNeeded( return isInCircumcircle(v, v2, v3, v4); } +template +TriInd Triangulation::edgeTriangle(const Edge edge) const +{ + TriInd iT = invalidIndex; + const TriInd start = m_vertTris[edge.v1()]; + TriInd currTri = start; + do + { + const Triangle& t = triangles[currTri]; + if(t.next(edge.v1()).second == edge.v2()) + { + iT = currTri; + break; + } + currTri = t.next(edge.v1()).first; + } while(currTri != start); + return iT; +} + +/* + * Contains a point in its diametral circle + * same as checking if the angle between v and edge end points is obtuse + */ +template +bool Triangulation::isEncroached( + const V2d& v, + const Edge edge) const +{ + const V2d v1 = + V2d::make(vertices[edge.v1()].x - v.x, vertices[edge.v1()].y - v.y); + const V2d v2 = + V2d::make(vertices[edge.v2()].x - v.x, vertices[edge.v2()].y - v.y); + return (v1.x * v2.x + v1.y * v2.y) < T(0); +} + +template +bool Triangulation::isBadTriangle( + const Triangle& tri, + const RefineTriangles::Enum refinement, + const T threshold) const +{ + const V2d& v1 = vertices[tri.vertices[0]]; + const V2d& v2 = vertices[tri.vertices[1]]; + const V2d& v3 = vertices[tri.vertices[2]]; + + const T twiceArea = std::abs(orient2D(v1, v2, v3)); + bool ans = false; + switch(refinement) + { + case RefineTriangles::ByAngle: { + T opoLenV1 = distance(v2, v3); + T opoLenV2 = distance(v1, v3); + T opoLenV3 = distance(v1, v2); + // Let opoLenV1 is the smallest edge length + if((opoLenV2 < opoLenV1) && (opoLenV2 < opoLenV3)) + { + std::swap(opoLenV1, opoLenV2); + } + else if(opoLenV3 < opoLenV1) + { + std::swap(opoLenV1, opoLenV3); + } + assert(opoLenV1 <= opoLenV3); + assert(opoLenV1 <= opoLenV2); + T samllestAngle = twiceArea / opoLenV3 / opoLenV2; + ans = samllestAngle < sin(threshold); + break; + } + case RefineTriangles::ByArea: { + T area = 0.5 * twiceArea; + ans = area > threshold; + break; + } + } + return ans; +} + +template +V2d +Triangulation::circumcenter(const Triangle& tri) const +{ + V2d a1 = vertices[tri.vertices[0]]; + V2d b1 = vertices[tri.vertices[1]]; + V2d a = vertices[tri.vertices[0]]; + V2d b = vertices[tri.vertices[1]]; + const V2d& c = vertices[tri.vertices[2]]; + const T denom = 0.5 / orient2D(c, a, b); + a.x -= c.x; + a.y -= c.y; + b.x -= c.x; + b.y -= c.y; + T oX = + c.x + + (b.y * (a.x * a.x + a.y * a.y) - a.y * (b.x * b.x + b.y * b.y)) * denom; + T oY = + c.y + + (a.x * (b.x * b.x + b.y * b.y) - b.x * (a.x * a.x + a.y * a.y)) * denom; + V2d v = V2d::make(oX, oY); + return v; +} + +/// Search in all fixed edges to find encroached edges, each fixed edge is +/// checked against its opposite vertices +/// Returns queue of encroached edges +template +EdgeQue Triangulation::detectEncroachedEdges() +{ + EdgeQue encroachedEdges; + for(EdgeUSet::const_iterator cit = fixedEdges.begin(); + cit != fixedEdges.end(); + ++cit) + { + const Edge edge = *cit; + const TriInd iT = edgeTriangle(edge); + const TriInd iTopo = edgeNeighbor(triangles[iT], edge.v1(), edge.v2()); + const Triangle& t = triangles[iT]; + const Triangle& tOpo = triangles[iTopo]; + VertInd v1 = opposedVertex(t, iTopo); + VertInd v2 = opposedVertex(tOpo, iT); + if(isEncroached(vertices[v1], edge) || isEncroached(vertices[v2], edge)) + { + encroachedEdges.push(edge); + } + } + return encroachedEdges; +} + +/// Search in all fixed edges to find encroached edges, each fixed edge is +/// checked against its opposite vertices and vertex v +/// Returns queue of encroached edges +template +EdgeQue +Triangulation::detectEncroachedEdges(const V2d& v) +{ + EdgeQue encroachedEdges; + for(EdgeUSet::const_iterator cit = fixedEdges.begin(); + cit != fixedEdges.end(); + ++cit) + { + const Edge edge = *cit; + if(isEncroached(v, edge)) + { + encroachedEdges.push(edge); + } + } + return encroachedEdges; +} + +template +std::pair +Triangulation::resolveEncroachedEdges( + EdgeQue encroachedEdges, + const V2d& v, + const bool validV, + const bool fillBadTriangles, + const RefineTriangles::Enum refinementConstrain, + const T badTriangleThreshold) +{ + IndexSizeType numOfSplits = 0; + std::vector badTriangles; + + while(!encroachedEdges.empty() && numOfSteinerPoints < maxSteinerPoints) + { + Edge edge = encroachedEdges.front(); + encroachedEdges.pop(); + if(!hasEdge(edge.v1(), edge.v2())) + { + continue; + } + TriInd iT = edgeTriangle(edge); + const Triangle& t = triangles[iT]; + VertInd i = splitEncroachedEdge( + edge, iT, edgeNeighbor(triangles[iT], edge.v1(), edge.v2())); + ++numOfSplits; + + TriInd start = m_vertTris[i]; + TriInd currTri = start; + do + { + const Triangle& t = triangles[currTri]; + if(fillBadTriangles && + isBadTriangle(t, refinementConstrain, badTriangleThreshold)) + { + badTriangles.push_back(currTri); + } + for(int i = 0; i < 3; ++i) + { + const Edge edge(t.vertices[i], t.vertices[cw(i)]); + if(fixedEdges.find(edge) == fixedEdges.end()) + continue; + const TriInd iT = currTri; + const TriInd iTopo = + edgeNeighbor(triangles[iT], edge.v1(), edge.v2()); + const Triangle& tOpo = triangles[iTopo]; + VertInd v1 = opposedVertex(t, iTopo); + VertInd v2 = opposedVertex(tOpo, iT); + if(isEncroached(vertices[v1], edge) || + isEncroached(vertices[v2], edge) || + (validV && isEncroached(v, edge))) + { + encroachedEdges.push(edge); + } + } + currTri = t.next(i).first; + } while(currTri != start); + } + return std::make_pair(badTriangles, numOfSplits); +} + +template +VertInd Triangulation::splitEncroachedEdge( + const Edge splitEdge, + const TriInd iT, + const TriInd iTopo) +{ + const VertInd iMid = static_cast(vertices.size()); + const V2d& start = vertices[splitEdge.v1()]; + const V2d& end = vertices[splitEdge.v2()]; + T split = T(0.5); + if(isSteinerVertex[splitEdge.v1()] || isSteinerVertex[splitEdge.v1()]) + { + // In Ruppert's paper, he used D(0.01) factor to divide edge length, but + // that introduces FP rounding erros, so it's avoided. + const T len = distance(start.x, start.y, end.x, end.y); + const T d = T(0.5) * len; + // Find the splitting distance. + T nearestPowerOfTwo = T(1); + while(d > nearestPowerOfTwo) + { + nearestPowerOfTwo *= T(2); + } + while(d < T(0.75) * nearestPowerOfTwo) + { + nearestPowerOfTwo *= T(0.5); + } + assert(abs(nearestPowerOfTwo - pow(2, round(log(d) / log(2.0)))) < 1e6); + split = nearestPowerOfTwo / len; + if(isSteinerVertex[splitEdge.v1()]) + split = T(1) - split; + } + V2d mid = V2d::make( + detail::lerp(start.x, end.x, split), + detail::lerp(start.y, end.y, split)); + + // split constraint edge that already exists in triangulation + if(fixedEdges.find(splitEdge) != fixedEdges.end()) + { + const Edge half1(splitEdge.v1(), iMid); + const Edge half2(iMid, splitEdge.v2()); + 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); + } + addNewVertex(mid, noNeighbor, true); + std::stack triStack = insertVertexOnEdge(iMid, iT, iTopo); + tryAddVertexToLocator(iMid); + ensureDelaunayByEdgeFlips(mid, iMid, triStack); + ++numOfSteinerPoints; + return iMid; +} + /* Flip edge between T and Topo: * * v4 | - old edge diff --git a/visualizer/main.cpp b/visualizer/main.cpp index fe57d922..fe74a0b9 100644 --- a/visualizer/main.cpp +++ b/visualizer/main.cpp @@ -192,6 +192,7 @@ public slots: CDT::VertexInsertionOrder::Auto, CDT::IntersectingConstraintEdges::Resolve, 1e-3); + m_cdt.maxSteinerPoints = 1000U; if(!m_points.empty()) { std::vector pts = @@ -219,12 +220,21 @@ public slots: else m_cdt.insertEdges(edges); } + if(m_isRemoveOuterAndHoles) + { + m_cdt.refineTriangles( + CDT::RefineTriangles::ByAngle, 20 / 180.0 * M_PI); m_cdt.eraseOuterTrianglesAndHoles(); + } else if(m_isRemoveOuter) m_cdt.eraseOuterTriangles(); else if(m_isHideSuperTri) + { + m_cdt.refineTriangles( + CDT::RefineTriangles::ByAngle, 20 / 180.0 * M_PI); m_cdt.eraseSuperTriangle(); + } const CDT::unordered_map tmp = CDT::EdgeToPiecesMapping(m_cdt.pieceToOriginals); const CDT::unordered_map >