Integrate Islam0mar's post-mirror Ruppert refinement fixes

Ports the two genuinely new commits from Islam0mar/CDT's
ruppert-delaunay-refinement branch that weren't yet in this mirror
(the rest of that branch's history was master commits pulled in via
periodic merges, already present here):

- aaac2dd "Refine after faces deletion + split fix"
- 3bda312 "avoid equidistant triangles in bad triangles"

Adapted by hand rather than cherry-picked, since this branch's own
subsequent refactoring (queue/criterion renames, isFlipNeeded/addNewVertex
signature changes on master) left little text in common with the
original diffs:

- isEncroachingOnEdge/isRefinementNeeded: use strict inequalities so a
  split landing exactly on a diametral circle/threshold boundary isn't
  re-flagged as encroached/bad, which could otherwise loop.
- splitEncroachedEdge: only apply the concentric-shell (power-of-two)
  split distance when the subsegment's Steiner endpoint is shared with
  another fixed edge at a small angle, instead of for every Steiner
  endpoint; drop an assert that silently truncated its operands to int.
- refineTriangles/resolveEncroachedEdges/splitEncroachedEdge take an
  optional toEraseOrNull set (matching the existing circumcenterOrNull
  pointer-optional convention in this file): triangles already marked
  for erasure are skipped as refinement candidates, and triangles that
  replace them via edge splits or vertex insertion inherit the mark.
  New collectSuperTriangle/collectOuterTriangles/
  collectOuterTrianglesAndHoles methods (paired with the now-public
  finalizeTriangulation) let a caller compute the erase set upfront and
  pass it to refineTriangles, so refinement doesn't waste the Steiner
  point budget on triangles about to be discarded (holes/outer
  triangles/super-triangle).
- Visualizer: collects the erase set for the selected finalize option
  before refining, and finalizes with the (possibly updated) set
  afterwards.

Verified: full CDT-tests suite (4787 assertions) and the CDT_visualizer
Qt build both pass; added ad hoc smoke tests exercising refineTriangles
standalone and combined with collectOuterTrianglesAndHoles + a hole,
confirming valid topology and no vertices placed outside the erased
region.

Author: artem-ogre

CDT/include/CDTUtils.hpp

@@ -278,7 +278,7 @@ bool isEncroachingOnEdge(
      * the angle between v and edge end points is obtuse
      */
     return (edgeStart.x - v.x) * (edgeEnd.x - v.x) +
-               (edgeStart.y - v.y) * (edgeEnd.y - v.y) <=
+               (edgeStart.y - v.y) * (edgeEnd.y - v.y) <
            T(0);
 }
 

CDT/include/Triangulation.h

@@ -698,12 +698,20 @@ class CDT_EXPORT Triangulation
      * @param refinement_constrain refinement strategy that is used to identify
      * bad triangles
      * @param refinementThreshold threshold value for refinement
+     * @param toEraseOrNull if not null, triangles already present in this set
+     * are skipped as refinement candidates and any triangles that replace them
+     * as a result of splitting are added to the set. Lets a caller combine
+     * refinement with one of the `eraseXXX` methods without wasting the vertex
+     * budget refining triangles that will be discarded anyway: pass in the
+     * result of the matching `collectXXX` method and later hand the (now
+     * updated) set to `finalizeTriangulation`.
      */
     void refineTriangles(
         VertInd maxVerticesToInsert,
         RefinementCriterion::Enum refinementCriterion =
             RefinementCriterion::SmallestAngle,
-        T refinementThreshold = 20 / 180.0 * M_PI);
+        T refinementThreshold = 20 / 180.0 * M_PI,
+        TriIndUSet* toEraseOrNull = NULL);
     /**
      * Erase triangles adjacent to super triangle
      *
@@ -719,6 +727,20 @@ class CDT_EXPORT Triangulation
      * @note supports overlapping or touching boundaries
      */
     void eraseOuterTrianglesAndHoles();
+    /**
+     * Collect triangles adjacent to super-triangle: same triangles that
+     * `eraseSuperTriangle` would remove.
+     * @note returns an empty set if custom geometry is used
+     */
+    TriIndUSet collectSuperTriangle() const;
+    /// Collect triangles outside of constrained boundary: same triangles that
+    /// `eraseOuterTriangles` would remove.
+    TriIndUSet collectOuterTriangles() const;
+    /**
+     * Collect triangles outside of constrained boundary and auto-detected
+     * holes: same triangles that `eraseOuterTrianglesAndHoles` would remove.
+     */
+    TriIndUSet collectOuterTrianglesAndHoles() const;
     /**
      * Call this method after directly setting custom super-geometry via
      * vertices and triangles members
@@ -802,6 +824,15 @@ class CDT_EXPORT Triangulation
     const TriIndVec& VertTrisInternal() const;
     /// @}
 
+    /**
+     * Remove super-triangle (if used) and triangles with specified indices.
+     * Adjust internal triangulation state accordingly.
+     * @param removedTriangles indices of triangles to remove
+     * @note pair with one of the `collectXXX` methods to combine erasing with
+     * `refineTriangles`
+     */
+    void finalizeTriangulation(const TriIndUSet& removedTriangles);
+
 private:
     /*____ Detail __*/
     void addSuperTriangle(const Box2d<T>& box);
@@ -940,8 +971,12 @@ class CDT_EXPORT Triangulation
         const V2d<T>* circumcenterOrNull = NULL,
         RefinementCriterion::Enum refinementCriterion =
             RefinementCriterion::SmallestAngle,
-        T badTriangleThreshold = T(0));
-    VertInd splitEncroachedEdge(Edge edge, VertInd steinerVerticesOffset);
+        T badTriangleThreshold = T(0),
+        TriIndUSet* toEraseOrNull = NULL);
+    VertInd splitEncroachedEdge(
+        Edge edge,
+        VertInd steinerVerticesOffset,
+        TriIndUSet* toEraseOrNull = NULL);
     void changeNeighbor(TriInd iT, TriInd oldNeighbor, TriInd newNeighbor);
     void changeNeighbor(
         TriInd iT,
@@ -966,12 +1001,6 @@ class CDT_EXPORT Triangulation
         IndexSizeType iB) const;
     TriInd addTriangle(const Triangle& t);
     TriInd addTriangle();
-    /**
-     * Remove super-triangle (if used) and triangles with specified indices.
-     * Adjust internal triangulation state accordingly.
-     * @removedTriangles indices of triangles to remove
-     */
-    void finalizeTriangulation(const TriIndUSet& removedTriangles);
     TriIndUSet growToBoundary(std::stack<TriInd> seeds) const;
     void fixEdge(const Edge& edge);
     void fixEdge(const Edge& edge, const Edge& originalEdge);

CDT/include/Triangulation.hpp

@@ -103,27 +103,47 @@ void Triangulation<T, TNearPointLocator>::eraseSuperTriangle()
 {
     if(m_superGeomType != SuperGeometryType::SuperTriangle)
         return;
+    finalizeTriangulation(collectSuperTriangle());
+}
+
+template <typename T, typename TNearPointLocator>
+void Triangulation<T, TNearPointLocator>::eraseOuterTriangles()
+{
+    finalizeTriangulation(collectOuterTriangles());
+}
+
+template <typename T, typename TNearPointLocator>
+void Triangulation<T, TNearPointLocator>::eraseOuterTrianglesAndHoles()
+{
+    finalizeTriangulation(collectOuterTrianglesAndHoles());
+}
+
+template <typename T, typename TNearPointLocator>
+TriIndUSet Triangulation<T, TNearPointLocator>::collectSuperTriangle() const
+{
     // find triangles adjacent to super-triangle's vertices
     TriIndUSet toErase;
+    if(m_superGeomType != SuperGeometryType::SuperTriangle)
+        return toErase;
     for(TriInd iT(0); iT < TriInd(triangles.size()); ++iT)
     {
         if(touchesSuperTriangle(triangles[iT]))
             toErase.insert(iT);
     }
-    finalizeTriangulation(toErase);
+    return toErase;
 }
 
 template <typename T, typename TNearPointLocator>
-void Triangulation<T, TNearPointLocator>::eraseOuterTriangles()
+TriIndUSet Triangulation<T, TNearPointLocator>::collectOuterTriangles() const
 {
     assert(m_vertTris[0] != noNeighbor);
     const std::stack<TriInd> seed(std::deque<TriInd>(1, m_vertTris[0]));
-    const TriIndUSet toErase = growToBoundary(seed);
-    finalizeTriangulation(toErase);
+    return growToBoundary(seed);
 }
 
 template <typename T, typename TNearPointLocator>
-void Triangulation<T, TNearPointLocator>::eraseOuterTrianglesAndHoles()
+TriIndUSet
+Triangulation<T, TNearPointLocator>::collectOuterTrianglesAndHoles() const
 {
     const std::vector<LayerDepth> triDepths = calculateTriangleDepths();
     TriIndUSet toErase;
@@ -133,7 +153,7 @@ void Triangulation<T, TNearPointLocator>::eraseOuterTrianglesAndHoles()
         if(triDepths[iT] % 2 == 0)
             toErase.insert(static_cast<TriInd>(iT));
     }
-    finalizeTriangulation(toErase);
+    return toErase;
 }
 
 /// Remap removing super-triangle: subtract 3 from vertices
@@ -1403,9 +1423,9 @@ bool Triangulation<T, TNearPointLocator>::isRefinementNeeded(
     switch(refinementCriterion)
     {
     case RefinementCriterion::SmallestAngle:
-        return smallestAngle(a, b, c) <= refinementThreshold;
+        return smallestAngle(a, b, c) < refinementThreshold;
     case RefinementCriterion::LargestArea:
-        return area(a, b, c) >= refinementThreshold;
+        return area(a, b, c) > refinementThreshold;
     }
     assert(false); // unreachable code
     return false;
@@ -1477,7 +1497,8 @@ TriIndVec Triangulation<T, TNearPointLocator>::resolveEncroachedEdges(
     const VertInd steinerVerticesOffset,
     const V2d<T>* const circumcenterOrNull,
     const RefinementCriterion::Enum refinementCriterion,
-    const T badTriangleThreshold)
+    const T badTriangleThreshold,
+    TriIndUSet* const toEraseOrNull)
 {
     std::vector<TriInd> badTriangles;
 
@@ -1490,16 +1511,18 @@ TriIndVec Triangulation<T, TNearPointLocator>::resolveEncroachedEdges(
             continue;
         }
         // split encroached edge
-        const VertInd iSplitVert =
-            splitEncroachedEdge(edge, steinerVerticesOffset);
+        const VertInd iSplitVert = splitEncroachedEdge(
+            edge, steinerVerticesOffset, toEraseOrNull);
         --remainingVertexBudget;
 
         const TriInd start = m_vertTris[iSplitVert];
         TriInd iT = start;
         do
         {
             const Triangle& t = triangles[iT];
-            if(circumcenterOrNull &&
+            const bool isMarkedForErasure =
+                toEraseOrNull && toEraseOrNull->count(iT);
+            if(circumcenterOrNull && !isMarkedForErasure &&
                isRefinementNeeded(t, refinementCriterion, badTriangleThreshold))
             {
                 badTriangles.push_back(iT);
@@ -1525,13 +1548,32 @@ TriIndVec Triangulation<T, TNearPointLocator>::resolveEncroachedEdges(
 template <typename T, typename TNearPointLocator>
 VertInd Triangulation<T, TNearPointLocator>::splitEncroachedEdge(
     const Edge edge,
-    const VertInd steinerVerticesOffset)
+    const VertInd steinerVerticesOffset,
+    TriIndUSet* const toEraseOrNull)
 {
     const V2d<T>& start = vertices[edge.v1()];
     const V2d<T>& end = vertices[edge.v2()];
+
+    TriInd iT, iTopo;
+    std::tie(iT, iTopo) = edgeTriangles(edge.v1(), edge.v2());
+    assert(iT != invalidIndex && iTopo != invalidIndex);
+
     T split = T(0.5);
-    // check if any of the split edge vertices are Steiner vertices
-    if(edge.v1() >= steinerVerticesOffset || edge.v2() >= steinerVerticesOffset)
+    // Use the concentric-shell splitting rule only when the edge is a
+    // subsegment terminating at a Steiner vertex that is itself an endpoint
+    // of another fixed (sub)segment meeting it at a small angle: splitting
+    // exactly in half in that case can lead to non-termination as the two
+    // segments keep encroaching on each other's ever-shrinking halves.
+    const VertInd v3 = opposedVertex(triangles[iT], iTopo);
+    const VertInd v4 = opposedVertex(triangles[iTopo], iT);
+    if((edge.v1() < steinerVerticesOffset &&
+        edge.v2() >= steinerVerticesOffset &&
+        (fixedEdges.find(Edge(v3, edge.v1())) != fixedEdges.end() ||
+         fixedEdges.find(Edge(v4, edge.v1())) != fixedEdges.end())) ||
+       (edge.v2() < steinerVerticesOffset &&
+        edge.v1() >= steinerVerticesOffset &&
+        (fixedEdges.find(Edge(v3, edge.v2())) != fixedEdges.end() ||
+         fixedEdges.find(Edge(v4, edge.v2())) != fixedEdges.end())))
     {
         // In Ruppert's paper, he used D(0.01) factor to divide edge length, but
         // that introduces FP rounding errors, so it's avoided.
@@ -1547,7 +1589,6 @@ VertInd Triangulation<T, TNearPointLocator>::splitEncroachedEdge(
         {
             nearestPowerOfTwo *= T(0.5);
         }
-        assert(abs(nearestPowerOfTwo - pow(2, round(log(d) / log(2.0)))) < 1e6);
         split = nearestPowerOfTwo / len;
         if(edge.v1() >= steinerVerticesOffset)
             split = T(1) - split;
@@ -1556,15 +1597,21 @@ VertInd Triangulation<T, TNearPointLocator>::splitEncroachedEdge(
     const V2d<T> mid = V2d<T>(
         detail::lerp(start.x, end.x, split),
         detail::lerp(start.y, end.y, split));
-    TriInd iT, iTopo;
-    std::tie(iT, iTopo) = edgeTriangles(edge.v1(), edge.v2());
-    assert(iT != invalidIndex && iTopo != invalidIndex);
 
     const VertInd iMid = addSplitEdgeVertex(mid, iT, iTopo);
     if(fixedEdges.find(edge) != fixedEdges.end())
     {
         splitFixedEdge(edge, iMid);
     }
+    // splitting reuses iT/iTopo for two of the four resulting triangles and
+    // appends the other two: propagate erasure marks to the new triangles
+    if(toEraseOrNull)
+    {
+        if(toEraseOrNull->count(iT))
+            toEraseOrNull->insert(TriInd(triangles.size() - 2));
+        if(toEraseOrNull->count(iTopo))
+            toEraseOrNull->insert(TriInd(triangles.size() - 1));
+    }
     return iMid;
 }
 
@@ -2393,7 +2440,8 @@ template <typename T, typename TNearPointLocator>
 void Triangulation<T, TNearPointLocator>::refineTriangles(
     const VertInd maxVerticesToInsert,
     const RefinementCriterion::Enum refinementCriterion,
-    const T refinementThreshold)
+    const T refinementThreshold,
+    TriIndUSet* const toEraseOrNull)
 {
     if(isFinalized())
     {
@@ -2412,8 +2460,8 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
     {
         const Edge edge = encroachedEdges.front();
         encroachedEdges.pop();
-        const VertInd iSplitVert =
-            splitEncroachedEdge(edge, steinerVerticesOffset);
+        const VertInd iSplitVert = splitEncroachedEdge(
+            edge, steinerVerticesOffset, toEraseOrNull);
         // if resulting halves are encroached, add them to the queue
         const Edge half1(edge.v1(), iSplitVert);
         if(isEdgeEncroached(half1))
@@ -2436,6 +2484,7 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
     {
         const Triangle& t = triangles[iT];
         if(!touchesSuperTriangle(t) &&
+           !(toEraseOrNull && toEraseOrNull->count(iT)) &&
            isRefinementNeeded(t, refinementCriterion, refinementThreshold))
         {
             badTriangles.push(iT);
@@ -2446,6 +2495,8 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
     {
         const TriInd iT = badTriangles.front();
         badTriangles.pop();
+        if(toEraseOrNull && toEraseOrNull->count(iT))
+            continue;
         const Triangle& badT = triangles[iT];
         if(!isRefinementNeeded(badT, refinementCriterion, refinementThreshold))
         {
@@ -2468,7 +2519,8 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
             steinerVerticesOffset,
             &triCircumenter,
             refinementCriterion,
-            refinementThreshold);
+            refinementThreshold,
+            toEraseOrNull);
         if(!remainingVertexBudget)
             break;
         if(!badTris.empty())
@@ -2482,17 +2534,45 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
             continue;
         }
 
-        --remainingVertexBudget;
+        // locate the triangle(s) the circumcenter falls into before spending
+        // the vertex budget: if it lands in a triangle that will be erased
+        // anyway (e.g. a hole or outer triangle) skip it instead of adding a
+        // pointless Steiner point there
         const VertInd iVert = static_cast<VertInd>(vertices.size());
         addNewVertex(triCircumenter, noNeighbor);
-        insertVertex(iVert);
+        const VertInd walkStart =
+            m_nearPtLocator.nearPoint(triCircumenter, vertices);
+        const array<TriInd, 2> trisAt =
+            walkingSearchTrianglesAt(iVert, walkStart);
+        if(toEraseOrNull &&
+           (toEraseOrNull->count(trisAt[0]) || toEraseOrNull->count(trisAt[1])))
+        {
+            vertices.pop_back();
+            m_vertTris.pop_back();
+            continue;
+        }
+
+        --remainingVertexBudget;
+        // re-does the same walk as above internally, but reuses proven
+        // insertion logic (fixed-edge split handling, callbacks) rather than
+        // duplicating it here
+        insertVertex(iVert, walkStart);
+        tryAddVertexToLocator(iVert);
+        if(toEraseOrNull)
+        {
+            if(toEraseOrNull->count(trisAt[0]))
+                toEraseOrNull->insert(TriInd(triangles.size() - 2));
+            if(toEraseOrNull->count(trisAt[1]))
+                toEraseOrNull->insert(TriInd(triangles.size() - 1));
+        }
 
         TriInd start = m_vertTris[iVert];
         TriInd currTri = start;
         do
         {
             const Triangle& t = triangles[currTri];
-            if(isRefinementNeeded(t, refinementCriterion, refinementThreshold))
+            if(!(toEraseOrNull && toEraseOrNull->count(currTri)) &&
+               isRefinementNeeded(t, refinementCriterion, refinementThreshold))
             {
                 badTriangles.push(currTri);
             }