#38: Frechet distance

.github/workflows/swift.yml

@@ -15,6 +15,7 @@ jobs:
     runs-on: macos-latest
 
     steps:
+    - uses: swift-actions/setup-swift@v1
     - uses: actions/checkout@v3
     - name: Clean Package
       run: swift package clean

Package.swift

@@ -1,4 +1,4 @@
-// swift-tools-version:5.7
+// swift-tools-version:5.9
 
 import PackageDescription
 

Sources/GISTools/Algorithms/FrechetDistance.swift

@@ -0,0 +1,91 @@
+#if !os(Linux)
+import CoreLocation
+#endif
+import Foundation
+
+/// Distance functions for ```LineString.frechetDistance```.
+public enum FrechetDistanceFunction {
+
+    /// Use the eudlidean distance.
+    case euclidean
+    /// Use the Haversine formula to account for global curvature.
+    case haversine
+    /// Use a rhumb line.
+    case rhumbLine
+    /// Use a custom distance function.
+    case other((Coordinate3D, Coordinate3D) -> CLLocationDistance)
+
+    func distance(between first: Coordinate3D, and second: Coordinate3D) -> CLLocationDistance {
+        switch self {
+        case .euclidean:
+            sqrt(pow(first.longitude - second.longitude, 2.0) + pow(first.latitude - second.latitude, 2.0))
+        case .haversine:
+            first.distance(from: second)
+        case .rhumbLine:
+            first.rhumbDistance(from: second)
+        case let .other(distanceFuntion):
+            distanceFuntion(first, second)
+        }
+    }
+
+}
+
+extension LineString {
+
+    /// Fréchet  distance between to geometries.
+    ///
+    /// - Parameters:
+    ///    - from: The other geometry of equal type.
+    ///    - distanceFunction: The algorithm to use for distance calculations.
+    ///    - segmentLength: Adds coordinates to the lines for improved matching (in meters).
+    ///
+    /// - Returns: The frechet distance between the two geometries.
+    public func frechetDistance(
+        from other: LineString,
+        distanceFunction: FrechetDistanceFunction = .haversine,
+        segmentLength: CLLocationDistance? = nil)
+        -> Double
+    {
+        var firstLine = self
+        var secondLine = other.projected(to: projection)
+
+        if let segmentLength {
+            firstLine = firstLine.evenlyDivided(segmentLength: segmentLength)
+            secondLine = secondLine.evenlyDivided(segmentLength: segmentLength)
+        }
+
+        let p = firstLine.allCoordinates
+        let q = secondLine.allCoordinates
+
+        var ca: [Double] = Array(repeating: -1.0, count: p.count * q.count)
+
+        func index(_ pI: Int, _ qI: Int) -> Int {
+            (pI * p.count) + qI
+        }
+
+        for i in 0 ..< p.count {
+            for j in 0 ..< q.count {
+                let distance = distanceFunction.distance(between: p[i], and: q[j])
+
+                ca[index(i, j)] = if i > 0, j > 0 {
+                    max([ca[index(i - 1, j)], ca[index(i - 1, j - 1)], ca[index(i, j - 1)]].min() ?? -1.0, distance)
+                }
+                else if i > 0, j == 0 {
+                    max(ca[index(i - 1, 0)], distance)
+                }
+                else if i == 0, j > 0 {
+                    max(ca[index(0, j - 1)], distance)
+                }
+                else if i == 0, j == 0 {
+                    distance
+                }
+                else {
+                    Double.infinity
+                }
+            }
+        }
+
+        return ca[index(p.count - 1, q.count - 1)]
+    }
+
+}

Sources/GISTools/Turf/LineChunk.swift → Sources/GISTools/Algorithms/LineChunk.swift

@@ -32,6 +32,18 @@ extension LineString {
         return MultiLineString(lineStrings) ?? MultiLineString()
     }
 
+    /// Divides a *LineString* into evenly spaced segments of a specified length.
+    /// If the line is shorter than the segment length then the original line is returned.
+    ///
+    /// - Parameter segmentLength: How long to make each segment, in meters
+    public func evenlyDivided(segmentLength: CLLocationDistance) -> LineString {
+        guard self.isValid, segmentLength > 0.0 else {
+            return LineString()
+        }
+
+        return LineString(chunked(segmentLength: segmentLength).lineSegments) ?? self
+    }
+
 }
 
 extension MultiLineString {
@@ -44,6 +56,16 @@ extension MultiLineString {
         MultiLineString(lineStrings.flatMap({ $0.chunked(segmentLength: segmentLength).lineStrings })) ?? self
     }
 
+    /// Divides a *MultiLineString* into evenly spaced segments of a specified length.
+    /// If the line is shorter than the segment length then the original line is returned.
+    ///
+    /// - Parameter segmentLength: How long to make each segment, in meters
+    public func evenlyDivided(segmentLength: CLLocationDistance) -> MultiLineString {
+        MultiLineString(lineStrings.map({
+            LineString($0.chunked(segmentLength: segmentLength).lineSegments) ?? $0
+        })) ?? self
+    }
+
 }
 
 extension Feature {

Tests/GISToolsTests/Algorithms/FrechetDistanceTests.swift

@@ -0,0 +1,30 @@
+#if !os(Linux)
+import CoreLocation
+#endif
+@testable import GISTools
+import XCTest
+
+final class FrechetDistanceTests: XCTestCase {
+
+    func testFrechetDistance4326() throws {
+        let point = Point(Coordinate3D(latitude: 0.0, longitude: 0.0))
+        let lineArc1 = try XCTUnwrap(point.lineArc(radius: 5000.0, bearing1: 20.0, bearing2: 60.0))
+        let lineArc2 = try XCTUnwrap(point.lineArc(radius: 6000.0, bearing1: 20.0, bearing2: 60.0))
+
+        let distanceHaversine = lineArc1.frechetDistance(from: lineArc2, distanceFunction: .haversine)
+        let distanceRhumbLine = lineArc1.frechetDistance(from: lineArc2, distanceFunction: .rhumbLine)
+
+        XCTAssertEqual(distanceHaversine, 1000.0, accuracy: 0.0001)
+        XCTAssertEqual(distanceRhumbLine, 1000.0, accuracy: 0.0001)
+    }
+
+    func testFrechetDistance3857() throws {
+        let point = Point(Coordinate3D(latitude: 0.0, longitude: 0.0)).projected(to: .epsg3857)
+        let lineArc1 = try XCTUnwrap(point.lineArc(radius: 5000.0, bearing1: 20.0, bearing2: 60.0))
+        let lineArc2 = try XCTUnwrap(point.lineArc(radius: 6000.0, bearing1: 20.0, bearing2: 60.0))
+
+        let distanceEucliden = lineArc1.frechetDistance(from: lineArc2, distanceFunction: .euclidean)
+        XCTAssertEqual(distanceEucliden, 1000.0, accuracy: 2.0)
+    }
+
+}

Tests/GISToolsTests/Turf/LineChunkTests.swift → Tests/GISToolsTests/Algorithms/LineChunkTests.swift

@@ -39,4 +39,19 @@ final class LineChunkTests: XCTestCase {
         XCTAssertEqual(chunks[0], lineString)
     }
 
+    func testEvenlyDivided() {
+        let a = Coordinate3D.zero
+        let b = a.destination(distance: 100.0, bearing: 90.0)
+        let line = LineString(unchecked: [a, b])
+        let dividedLine = line.evenlyDivided(segmentLength: 1.0)
+
+        XCTAssertEqual(line.allCoordinates.count, 2)
+        XCTAssertEqual(dividedLine.allCoordinates.count, 101)
+
+        for (first, second) in dividedLine.allCoordinates.overlappingPairs() {
+            guard let second else { break }
+            XCTAssertEqual(first.distance(from: second), 1.0, accuracy: 0.0001)
+        }
+    }
+
 }