#6182 Improve calculate_circle_center (#6192)

* add tests for calculate_circle_center - the first one succeeds with the current impl

* fix calculate_circle_center

* comment cleanup

* clippy

* comment format

* update circle_three_point sim test snapshot for slight floating point changes introduced by calculate_circle_center refactor

* use TAU instead of 2 * PI

* clippy

rust/kcl-lib/src/std/utils.rs

@@ -234,40 +234,37 @@ pub fn is_on_circumference(center: Point2d, point: Point2d, radius: f64) -> bool
     (distance_squared - radius.powi(2)).abs() < 1e-9
 }
 
-// Calculate the center of 3 points
+// Calculate the center of 3 points using an algebraic method
-// To calculate the center of the 3 point circle 2 perpendicular lines are created
+// Handles if 3 points lie on the same line (collinear) by returning the average of the points (could return None instead..)
-// These perpendicular lines will intersect at the center of the circle.
 pub fn calculate_circle_center(p1: [f64; 2], p2: [f64; 2], p3: [f64; 2]) -> [f64; 2] {
-    // y2 - y1
+    let (x1, y1) = (p1[0], p1[1]);
-    let y_2_1 = p2[1] - p1[1];
+    let (x2, y2) = (p2[0], p2[1]);
-    // y3 - y2
+    let (x3, y3) = (p3[0], p3[1]);
-    let y_3_2 = p3[1] - p2[1];
-    // x2 - x1
-    let x_2_1 = p2[0] - p1[0];
-    // x3 - x2
-    let x_3_2 = p3[0] - p2[0];
 
-    // Slope of two perpendicular lines
+    // Compute the determinant d = 2 * (x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
-    let slope_a = y_2_1 / x_2_1;
+    // Visually d is twice the area of the triangle formed by the points,
-    let slope_b = y_3_2 / x_3_2;
+    // also the same as: cross(p2 - p1, p3 - p1)
+    let d = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));
 
-    // Values for line intersection
+    // If d is nearly zero, the points are collinear, and a unique circle cannot be defined.
-    // y1 - y3
+    if d.abs() < f64::EPSILON {
-    let y_1_3 = p1[1] - p3[1];
+        return [(x1 + x2 + x3) / 3.0, (y1 + y2 + y3) / 3.0];
-    // x1 + x2
+    }
-    let x_1_2 = p1[0] + p2[0];
-    // x2 + x3
-    let x_2_3 = p2[0] + p3[0];
-    // y1 + y2
-    let y_1_2 = p1[1] + p2[1];
 
-    // Solve for the intersection of these two lines
+    // squared lengths
-    let numerator = (slope_a * slope_b * y_1_3) + (slope_b * x_1_2) - (slope_a * x_2_3);
+    let p1_sq = x1 * x1 + y1 * y1;
-    let x = numerator / (2.0 * (slope_b - slope_a));
+    let p2_sq = x2 * x2 + y2 * y2;
+    let p3_sq = x3 * x3 + y3 * y3;
 
-    let y = ((-1.0 / slope_a) * (x - (x_1_2 / 2.0))) + (y_1_2 / 2.0);
+    // This formula is derived from the circle equations:
-
+    //   (x - cx)^2 + (y - cy)^2 = r^2
-    [x, y]
+    // All 3 points will satisfy this equation, so we have 3 equations. Radius can be eliminated
+    // by subtracting one of the equations from the other two and the remaining 2 equations can
+    // be solved for cx and cy.
+    [
+        (p1_sq * (y2 - y3) + p2_sq * (y3 - y1) + p3_sq * (y1 - y2)) / d,
+        (p1_sq * (x3 - x2) + p2_sq * (x1 - x3) + p3_sq * (x2 - x1)) / d,
+    ]
 }
 
 pub struct CircleParams {
@@ -286,9 +283,11 @@ pub fn calculate_circle_from_3_points(points: [Point2d; 3]) -> CircleParams {
 #[cfg(test)]
 mod tests {
     // Here you can bring your functions into scope
+    use approx::assert_relative_eq;
     use pretty_assertions::assert_eq;
+    use std::f64::consts::TAU;
 
-    use super::{get_x_component, get_y_component, Angle};
+    use super::{calculate_circle_center, get_x_component, get_y_component, Angle};
     use crate::SourceRange;
 
     static EACH_QUAD: [(i32, [i32; 2]); 12] = [
@@ -453,6 +452,75 @@ mod tests {
         assert_eq!(angle_start.to_degrees().round(), 0.0);
         assert_eq!(angle_end.to_degrees().round(), 180.0);
     }
+
+    #[test]
+    fn test_calculate_circle_center() {
+        const EPS: f64 = 1e-4;
+
+        // Test: circle center = (4.1, 1.9)
+        let p1 = [1.0, 2.0];
+        let p2 = [4.0, 5.0];
+        let p3 = [7.0, 3.0];
+        let center = calculate_circle_center(p1, p2, p3);
+        assert_relative_eq!(center[0], 4.1, epsilon = EPS);
+        assert_relative_eq!(center[1], 1.9, epsilon = EPS);
+
+        // Tests: Generate a few circles and test its points
+        let center = [3.2, 0.7];
+        let radius_array = [0.001, 0.01, 0.6, 1.0, 5.0, 60.0, 500.0, 2000.0, 400_000.0];
+        let points_array = [[0.0, 0.33, 0.66], [0.0, 0.1, 0.2], [0.0, -0.1, 0.1], [0.0, 0.5, 0.7]];
+
+        let get_point = |radius: f64, t: f64| {
+            let angle = t * TAU;
+            [center[0] + radius * angle.cos(), center[1] + radius * angle.sin()]
+        };
+
+        for radius in radius_array {
+            for point in points_array {
+                let p1 = get_point(radius, point[0]);
+                let p2 = get_point(radius, point[1]);
+                let p3 = get_point(radius, point[2]);
+                let c = calculate_circle_center(p1, p2, p3);
+                assert_relative_eq!(c[0], center[0], epsilon = EPS);
+                assert_relative_eq!(c[1], center[1], epsilon = EPS);
+            }
+        }
+
+        // Test: Equilateral triangle
+        let p1 = [0.0, 0.0];
+        let p2 = [1.0, 0.0];
+        let p3 = [0.5, 3.0_f64.sqrt() / 2.0];
+        let center = calculate_circle_center(p1, p2, p3);
+        assert_relative_eq!(center[0], 0.5, epsilon = EPS);
+        assert_relative_eq!(center[1], 1.0 / (2.0 * 3.0_f64.sqrt()), epsilon = EPS);
+
+        // Test: Collinear points (should return the average of the points)
+        let p1 = [0.0, 0.0];
+        let p2 = [1.0, 0.0];
+        let p3 = [2.0, 0.0];
+        let center = calculate_circle_center(p1, p2, p3);
+        assert_relative_eq!(center[0], 1.0, epsilon = EPS);
+        assert_relative_eq!(center[1], 0.0, epsilon = EPS);
+
+        // Test: Points forming a circle with radius = 1
+        let p1 = [0.0, 0.0];
+        let p2 = [0.0, 2.0];
+        let p3 = [2.0, 0.0];
+        let center = calculate_circle_center(p1, p2, p3);
+        assert_relative_eq!(center[0], 1.0, epsilon = EPS);
+        assert_relative_eq!(center[1], 1.0, epsilon = EPS);
+
+        // Test: Integer coordinates
+        let p1 = [0.0, 0.0];
+        let p2 = [0.0, 6.0];
+        let p3 = [6.0, 0.0];
+        let center = calculate_circle_center(p1, p2, p3);
+        assert_relative_eq!(center[0], 3.0, epsilon = EPS);
+        assert_relative_eq!(center[1], 3.0, epsilon = EPS);
+        // Verify radius (should be 3 * sqrt(2))
+        let radius = ((center[0] - p1[0]).powi(2) + (center[1] - p1[1]).powi(2)).sqrt();
+        assert_relative_eq!(radius, 3.0 * 2.0_f64.sqrt(), epsilon = EPS);
+    }
 }
 
 pub type Coords2d = [f64; 2];

rust/kcl-lib/tests/circle_three_point/artifact_commands.snap

@@ -85,7 +85,7 @@ description: Artifact commands circle_three_point.kcl
       "path": "[uuid]",
       "to": {
         "x": 30.00594901040716,
-        "y": 19.749999999999996,
+        "y": 19.75,
         "z": 0.0
       }
     }
@@ -109,7 +109,7 @@ description: Artifact commands circle_three_point.kcl
           "x": 24.75,
           "y": 19.75
         },
-        "radius": 5.255949010407163,
+        "radius": 5.25594901040716,
         "start": {
           "unit": "degrees",
           "value": 0.0