Adding GeometriaAlgorithms target with 2D ellipse closest point appro…

…ximation algorithm

Package.swift

@@ -25,10 +25,22 @@ let geometriaTestsTarget: Target = .testTarget(
     dependencies: ["Geometria"],
     swiftSettings: []
 )
+let geometriaAlgorithmsTarget: Target = .target(
+    name: "GeometriaAlgorithms",
+    dependencies: geometriaDependencies + ["Geometria"],
+    swiftSettings: []
+)
+let geometriaAlgorithmsTestTarget: Target = .testTarget(
+    name: "GeometriaAlgorithmsTests",
+    dependencies: geometriaDependencies + ["GeometriaAlgorithms"],
+    swiftSettings: []
+)
 
 if ProcessInfo.processInfo.environment["REPORT_BUILD_TIME"] == "YES" {
     geometriaTarget.swiftSettings?.append(contentsOf: reportingSwiftSettings)
     geometriaTestsTarget.swiftSettings?.append(contentsOf: reportingSwiftSettings)
+    geometriaAlgorithmsTarget.swiftSettings?.append(contentsOf: reportingSwiftSettings)
+    geometriaAlgorithmsTestTarget.swiftSettings?.append(contentsOf: reportingSwiftSettings)
 }
 
 let package = Package(
@@ -44,5 +56,7 @@ let package = Package(
     targets: [
         geometriaTarget,
         geometriaTestsTarget,
+        geometriaAlgorithmsTarget,
+        geometriaAlgorithmsTestTarget,
     ]
 )

Sources/GeometriaAlgorithms/2D/Ellipse2.swift

@@ -0,0 +1,76 @@
+import Numerics
+import Geometria
+
+public extension Ellipse2 where Vector: VectorReal {
+    /// Returns the point on a given angle, in radians, on this ellipse.
+    @inlinable
+    func pointOnAngle(_ angleInRadians: Scalar) -> Vector {
+        let c = Scalar.cos(angleInRadians)
+        let s = Scalar.sin(angleInRadians)
+
+        return center + .init(x: c, y: s) * radius
+    }
+
+    /// Performs an approximation-type search to find the closest point on the
+    /// ellipse's outer surface to a given point.
+    ///
+    /// Each step of the search samples `samples` points on the ellipse, using
+    /// the closest point on the samples to perform a deeper search around the
+    /// sample, up to a depth of `maxDepth`, or until a mean tolerance of
+    /// `tolerance` between samples is achieved.
+    ///
+    /// - complexity: `O(samples * maxDepth)`, for the worst case of the search.
+    func approximateClosestPoint(
+        to vector: Vector,
+        tolerance: Scalar,
+        samples: Int,
+        maxDepth: Int
+    ) -> Vector where Scalar: Strideable, Scalar.Stride == Scalar {
+
+        precondition(samples > 0, "samples > 0")
+        precondition(maxDepth > 0, "maxDepth > 0")
+
+        var closestInRadians: (distanceSquared: Scalar, angle: Scalar) = (.infinity, 0)
+
+        let samplesAsScalar = Scalar(samples)
+
+        var searchStart: Scalar = 0
+        var searchEnd: Scalar = .pi * 2
+        var searchStep: Scalar = .pi * 2 / samplesAsScalar
+
+        outer:
+        for _ in 0..<maxDepth {
+            var lastPointDistance: Scalar = .infinity
+
+            for radians in stride(from: searchStart, to: searchEnd, by: searchStep) {
+                let point = pointOnAngle(radians)
+                let distanceSquared = point.distanceSquared(to: vector)
+
+                if distanceSquared < closestInRadians.distanceSquared {
+                    closestInRadians = (distanceSquared, radians)
+                }
+
+                // Check if we are increasing in distance, instead of decreasing,
+                // and quit early; unless we are in the first depth of search, in
+                // which case span the entire sample space, since leaving it out
+                // might result in the wrong quadrant being picked
+                if lastPointDistance < distanceSquared {
+                    break
+                }
+
+                lastPointDistance = distanceSquared
+            }
+
+            if searchStep < tolerance {
+                break
+            }
+
+            // Narrow search
+            searchStep /= samplesAsScalar
+            searchStart = closestInRadians.angle - searchStep * 10
+            searchEnd = closestInRadians.angle + searchStep * 10
+        }
+
+        return pointOnAngle(closestInRadians.angle)
+    }
+}

Tests/GeometriaAlgorithmsTests/2D/Ellipse2Tests.swift

@@ -0,0 +1,41 @@
+import XCTest
+import Geometria
+
+@testable import GeometriaAlgorithms
+
+class Ellipse2Tests: XCTestCase {
+    func testApproximateClosestPoint_centerPoint() {
+        let sut = makeSut(center: .zero, radius: .init(x: 10, y: 10))
+        let point = Vector2D(x: 0, y: 0)
+
+        let result = sut.approximateClosestPoint(to: point, tolerance: 0.0001, samples: 20, maxDepth: 4)
+
+        assertEqual(result, .init(x: 9.510565162951535, y: 3.090169943749474), accuracy: 0.001)
+    }
+
+    func testApproximateClosestPoint_circleEquivalent() {
+        let sut = makeSut(center: .zero, radius: .init(x: 10, y: 10))
+        let circle = Circle2D(center: .zero, radius: 10)
+        let point = Vector2D(x: 3, y: 4)
+        let onCircle = circle.project(point)
+
+        let result = sut.approximateClosestPoint(to: point, tolerance: 0.0001, samples: 20, maxDepth: 4)
+
+        assertEqual(result, onCircle, accuracy: 0.001)
+    }
+
+    func testApproximateClosestPoint_nonCircleEquivalent() {
+        let sut = makeSut(center: .init(x: 2, y: 3), radius: .init(x: 12, y: 7))
+        let point = Vector2D(x: 5, y: 7)
+
+        let result = sut.approximateClosestPoint(to: point, tolerance: 0.0001, samples: 20, maxDepth: 5)
+
+        assertEqual(result, .init(x: 5.476651310470384, y: 9.699778130449972), accuracy: 0.001)
+    }
+}
+
+// MARK: - Test internals
+
+private func makeSut(center: Vector2D, radius: Vector2D) -> Ellipse2D {
+    .init(center: center, radius: radius)
+}