Ruppert's Refinement Algorithm.

.github/workflows/ci_builds.yml

@@ -19,17 +19,17 @@ jobs:
             os: ubuntu-latest
           }
         - {
-            name: "macOS",
+            name: "macOS", 
             os: macos-latest
           }
 
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v3
       with:
         submodules: true
 
     - name: Set up Python
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v2
       with:
         python-version: '3.x'
 
@@ -39,7 +39,7 @@ jobs:
         pip install conan
         conan --version
         conan profile detect
-
+    
     - name: "Build: compiled, shared, boost"
       run: conan create CDT/ -o as_compiled_library=True -o shared=True -o use_boost=True -o enable_testing=True --build missing
     - name: "Build: compiled, shared, without boost"

CDT/conanfile.py

@@ -29,8 +29,8 @@ class CDTConan(ConanFile):
 
     def requirements(self):
         if self.options.use_boost:
-            self.requires("boost/1.83.0")
-        self.requires("catch2/3.4.0")
+            self.requires("boost/1.81.0")
+        self.requires("catch2/3.3.1")
 
     def configure(self):
         if self.options.use_boost:

CDT/extras/VerifyTopology.h

@@ -74,17 +74,6 @@ inline bool verifyTopology(const CDT::Triangulation<T, TNearPointLocator>& cdt)
     return true;
 }
 
-/// Check that each vertex has a neighbor triangle
-template <typename T, typename TNearPointLocator>
-inline bool eachVertexHasNeighborTriangle(
-    const CDT::Triangulation<T, TNearPointLocator>& cdt)
-{
-    for(const auto& vt : cdt.VertTrisInternal())
-        if(vt == noNeighbor)
-            return false;
-    return true;
-}
-
 } // namespace CDT
 
 #endif

CDT/include/CDTUtils.h

@@ -44,6 +44,7 @@ typedef char couldnt_parse_cxx_standard[-1]; ///< Error: couldn't parse standard
 #include <cassert>
 #include <cmath>
 #include <limits>
+#include <queue>
 #include <vector>
 
 #ifdef CDT_USE_STRONG_TYPING
@@ -243,6 +244,7 @@ inline Edge edge_make(VertInd iV1, VertInd iV2)
 }
 
 typedef std::vector<Edge> EdgeVec;                ///< Vector of edges
+typedef std::queue<Edge> EdgeQueue;                 ///< Queue of edges
 typedef unordered_set<Edge> EdgeUSet;             ///< Hash table of edges
 typedef unordered_set<TriInd> TriIndUSet;         ///< Hash table of triangles
 typedef unordered_map<TriInd, TriInd> TriIndUMap; ///< Triangle hash map
@@ -324,7 +326,6 @@ struct CDT_EXPORT PtTriLocation
         OnEdge1,
         OnEdge2,
         OnEdge3,
-        OnVertex,
     };
 };
 
@@ -418,6 +419,14 @@ CDT_EXPORT bool isInCircumcircle(
 CDT_EXPORT CDT_INLINE_IF_HEADER_ONLY bool
 verticesShareEdge(const TriIndVec& aTris, const TriIndVec& bTris);
 
+/// Vector's length
+template <typename T>
+CDT_EXPORT T length(const V2d<T>& v);
+
+/// Vector's squared length
+template <typename T>
+CDT_EXPORT T lengthSquared(const V2d<T>& v);
+
 /// Distance between two 2D points
 template <typename T>
 CDT_EXPORT T distance(const V2d<T>& a, const V2d<T>& b);
@@ -426,8 +435,24 @@ CDT_EXPORT T distance(const V2d<T>& a, const V2d<T>& b);
 template <typename T>
 CDT_EXPORT T distanceSquared(const V2d<T>& a, const V2d<T>& b);
 
-/// Check if any of triangle's vertices belongs to a super-triangle
-CDT_INLINE_IF_HEADER_ONLY bool touchesSuperTriangle(const Triangle& t);
+/// Check if vertex V is encroaching on diametral circle of an edge
+template <typename T>
+CDT_EXPORT bool isEncroachingOnEdge(
+    const V2d<T>& v,
+    const V2d<T>& edgeStart,
+    const V2d<T>& edgeEnd);
+
+/// Position of ABC triangle circumcenter
+template <typename T>
+CDT_EXPORT V2d<T> circumcenter(V2d<T> a, V2d<T> b, const V2d<T>& c);
+
+/// Doubled surface area of a triangle ABC
+template <typename T>
+CDT_EXPORT T doubledArea(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c);
+
+/// Surface area of a triangle ABC
+template <typename T>
+CDT_EXPORT T area(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c);
 
 } // namespace CDT
 

CDT/include/CDTUtils.hpp

@@ -84,14 +84,12 @@ CDT_INLINE_IF_HEADER_ONLY Index cw(Index i)
 
 CDT_INLINE_IF_HEADER_ONLY bool isOnEdge(const PtTriLocation::Enum location)
 {
-    return location == PtTriLocation::OnEdge1 ||
-           location == PtTriLocation::OnEdge2 ||
-           location == PtTriLocation::OnEdge3;
+    return location > PtTriLocation::Outside;
 }
 
 CDT_INLINE_IF_HEADER_ONLY Index edgeNeighbor(const PtTriLocation::Enum location)
 {
-    assert(isOnEdge(location));
+    assert(location >= PtTriLocation::OnEdge1);
     return static_cast<Index>(location - PtTriLocation::OnEdge1);
 }
 
@@ -140,18 +138,12 @@ PtTriLocation::Enum locatePointTriangle(
     if(edgeCheck == PtLineLocation::Right)
         return PtTriLocation::Outside;
     if(edgeCheck == PtLineLocation::OnLine)
-    {
-        result = (result == PtTriLocation::Inside) ? PtTriLocation::OnEdge2
-                                                   : PtTriLocation::OnVertex;
-    }
+        result = PtTriLocation::OnEdge2;
     edgeCheck = locatePointLine(p, v3, v1);
     if(edgeCheck == PtLineLocation::Right)
         return PtTriLocation::Outside;
     if(edgeCheck == PtLineLocation::OnLine)
-    {
-        result = (result == PtTriLocation::Inside) ? PtTriLocation::OnEdge3
-                                                   : PtTriLocation::OnVertex;
-    }
+        result = PtTriLocation::OnEdge3;
     return result;
 }
 
@@ -163,7 +155,6 @@ CDT_INLINE_IF_HEADER_ONLY Index opoNbr(const Index vertIndex)
         return Index(2);
     if(vertIndex == Index(2))
         return Index(0);
-    assert(false && "Invalid vertex index");
     throw std::runtime_error("Invalid vertex index");
 }
 
@@ -175,7 +166,6 @@ CDT_INLINE_IF_HEADER_ONLY Index opoVrt(const Index neighborIndex)
         return Index(0);
     if(neighborIndex == Index(2))
         return Index(1);
-    assert(false && "Invalid neighbor index");
     throw std::runtime_error("Invalid neighbor index");
 }
 
@@ -285,34 +275,72 @@ bool verticesShareEdge(const TriIndVec& aTris, const TriIndVec& bTris)
 }
 
 template <typename T>
-T distanceSquared(const T ax, const T ay, const T bx, const T by)
+T lengthSquared(const T x, const T y)
 {
-    const T dx = bx - ax;
-    const T dy = by - ay;
-    return dx * dx + dy * dy;
+    return x * x + y * y;
 }
 
 template <typename T>
-T distance(const T ax, const T ay, const T bx, const T by)
+T lengthSquared(const V2d<T>& v)
 {
-    return std::sqrt(distanceSquared(ax, ay, bx, by));
+    return lengthSquared(v.x, v.y);
 }
 
 template <typename T>
-T distance(const V2d<T>& a, const V2d<T>& b)
+T length(const V2d<T>& v)
 {
-    return distance(a.x, a.y, b.x, b.y);
+    return std::sqrt(lengthSquared(v));
 }
 
 template <typename T>
 T distanceSquared(const V2d<T>& a, const V2d<T>& b)
 {
-    return distanceSquared(a.x, a.y, b.x, b.y);
+    return lengthSquared(b.x - a.x, b.y - a.y);
+}
+
+template <typename T>
+T distance(const V2d<T>& a, const V2d<T>& b)
+{
+    return std::sqrt(distanceSquared(a, b));
+}
+
+template <typename T>
+bool isEncroachingOnEdge(
+    const V2d<T>& v,
+    const V2d<T>& edgeStart,
+    const V2d<T>& edgeEnd)
+{
+    /*
+     * Contains a point in its diametral circle:
+     * 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) <
+           T(0);
+}
+
+template <typename T>
+V2d<T> circumcenter(V2d<T> a, V2d<T> b, const V2d<T>& c)
+{
+    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;
+    const T aLenSq = lengthSquared(a), bLenSq = lengthSquared(b);
+    return V2d<T>::make(
+        c.x + (b.y * aLenSq - a.y * bLenSq) * denom,
+        c.y + (a.x * bLenSq - b.x * aLenSq) * denom);
+}
+
+template <typename T>
+T doubledArea(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c)
+{
+    return std::abs(orient2D(a, b, c));
 }
 
-bool touchesSuperTriangle(const Triangle& t)
+template <typename T>
+T area(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c)
 {
-    return t.vertices[0] < 3 || t.vertices[1] < 3 || t.vertices[2] < 3;
+    return doubledArea(a, b, c) / T(2);
 }
 
 } // namespace CDT