From 9740eee4b173d5e0168216a0c00c8f1ea4d8f321 Mon Sep 17 00:00:00 2001
From: artem-ogre <artem.ogre@gmail.com>
Date: Mon, 21 Aug 2023 11:40:44 +0200
Subject: [PATCH] #7 Refactor calculating triangle circumcenter

---
 CDT/include/CDTUtils.h        | 12 ++++++++++++
 CDT/include/CDTUtils.hpp      | 34 ++++++++++++++++++++++++---------
 CDT/include/Triangulation.h   |  1 -
 CDT/include/Triangulation.hpp | 36 +++++++++--------------------------
 4 files changed, 46 insertions(+), 37 deletions(-)

diff --git a/CDT/include/CDTUtils.h b/CDT/include/CDTUtils.h
index 3452f46..dda19c4 100644
--- a/CDT/include/CDTUtils.h
+++ b/CDT/include/CDTUtils.h
@@ -419,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);
@@ -434,6 +442,10 @@ CDT_EXPORT bool isEncroachingOnEdge(
     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);
+
 } // namespace CDT
 
 #ifndef CDT_USE_AS_COMPILED_LIBRARY
diff --git a/CDT/include/CDTUtils.hpp b/CDT/include/CDTUtils.hpp
index f068cac..62605ad 100644
--- a/CDT/include/CDTUtils.hpp
+++ b/CDT/include/CDTUtils.hpp
@@ -275,29 +275,33 @@ 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>
@@ -315,4 +319,16 @@ bool isEncroachingOnEdge(
            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);
+}
+
 } // namespace CDT
diff --git a/CDT/include/Triangulation.h b/CDT/include/Triangulation.h
index 3911abe..119ae6c 100644
--- a/CDT/include/Triangulation.h
+++ b/CDT/include/Triangulation.h
@@ -525,7 +525,6 @@ class CDT_EXPORT Triangulation
         const Triangle& tri,
         RefineTriangles::Enum refinement,
         T threshold) const;
-    V2d<T> 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
diff --git a/CDT/include/Triangulation.hpp b/CDT/include/Triangulation.hpp
index 3f08e30..9b05d46 100644
--- a/CDT/include/Triangulation.hpp
+++ b/CDT/include/Triangulation.hpp
@@ -1327,30 +1327,6 @@ bool Triangulation<T, TNearPointLocator>::isBadTriangle(
     return ans;
 }
 
-template <typename T, typename TNearPointLocator>
-V2d<T>
-Triangulation<T, TNearPointLocator>::circumcenter(const Triangle& tri) const
-{
-    V2d<T> a1 = vertices[tri.vertices[0]];
-    V2d<T> b1 = vertices[tri.vertices[1]];
-    V2d<T> a = vertices[tri.vertices[0]];
-    V2d<T> b = vertices[tri.vertices[1]];
-    const V2d<T>& 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<T> v = V2d<T>::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
@@ -1478,7 +1454,7 @@ VertInd Triangulation<T, TNearPointLocator>::splitEncroachedEdge(
     {
         // 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 len = distance(start, end);
         const T d = T(0.5) * len;
         // Find the splitting distance.
         T nearestPowerOfTwo = T(1);
@@ -2341,7 +2317,10 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
 
         if(isBadTriangle(t, refinementConstrain, threshold))
         {
-            const V2d<T> vert = circumcenter(t);
+            const V2d<T> vert = circumcenter(
+                vertices[t.vertices[0]],
+                vertices[t.vertices[1]],
+                vertices[t.vertices[2]]);
             if(locatePointTriangle(
                    vert, vertices[0], vertices[1], vertices[2]) !=
                PtTriLocation::Outside)
@@ -2361,7 +2340,10 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
         {
             continue;
         }
-        const V2d<T> vert = circumcenter(t);
+        const V2d<T> vert = circumcenter(
+            vertices[t.vertices[0]],
+            vertices[t.vertices[1]],
+            vertices[t.vertices[2]]);
         if(locatePointTriangle(vert, vertices[0], vertices[1], vertices[2]) ==
            PtTriLocation::Outside)
         {