Discrete Frechet Distance (#242)

* discrete frechet distance using dynamic programming

* generalized the types + some cleaning up

hgeometry/hgeometry.cabal

@@ -367,6 +367,9 @@ library
         HGeometry.PolyLine.Simplification.DouglasPeucker
         HGeometry.PolyLine.Simplification.ImaiIri
 
+        HGeometry.PolyLine.Frechet.Discrete
+
+
 --
         HGeometry.ConvexHull
         HGeometry.ConvexHull.GrahamScan
@@ -705,6 +708,7 @@ test-suite hspec
       Polygon.Triangulation.MakeMonotoneSpec
       Polygon.VerticesSpec
       PolyLine.Simplification.ImaiIriSpec
+      PolyLine.Frechet.DiscreteSpec
       SegmentTreeSpec
       SegmentTree.R2Spec
       RangeTreeSpec

hgeometry/src/HGeometry/PolyLine/Frechet/Discrete.hs

@@ -0,0 +1,60 @@
+--------------------------------------------------------------------------------
+-- |
+-- Module      :  HGeometry.PolyLine.Frechet.Discrete
+-- Copyright   :  (C) Frank Staals
+-- License     :  see the LICENSE file
+-- Maintainer  :  Frank Staals
+--
+-- Computes the discrete frechet distance
+--
+--------------------------------------------------------------------------------
+module HGeometry.PolyLine.Frechet.Discrete
+  ( frechetDistanceWith
+  ) where
+
+import Control.Lens
+import Data.Array
+import Hiraffe.Graph
+
+--------------------------------------------------------------------------------
+
+-- | Computes the discrete frechet distance with respect to the given distance function.
+--
+-- O(nm*(I+T)), where I is the time to index one of the vertices and T is the time to
+-- evaluate the distance function.
+frechetDistanceWith :: ( HasVertices' polyLine
+                       , HasVertices' polyLine'
+                       , Ord r
+                       , VertexIx polyLine ~ Int, VertexIx polyLine' ~ Int
+                       )
+                    => (Vertex polyLine -> Vertex polyLine' -> r)
+                    -> polyLine -> polyLine'
+                    -> r
+frechetDistanceWith dist' polyP polyQ = fd (n,m)
+  where
+    n = numVertices polyP - 1
+    m = numVertices polyQ - 1
+    dist i j = dist' (polyP^?!vertexAt i) (polyQ^?!vertexAt j)
+
+    -- use dynamic programming
+    fd = memo ((0, 0), (n, m)) fd'
+
+    fd' (0,0) = dist 0 0
+    fd' (0,q) = dist 0 q `max` fd (0, q-1)
+    fd' (p,0) = dist p 0 `max` fd (p-1, 0)
+    fd' (p,q) = dist p q `max` (fd (p, q-1) `min` fd (p-1, q) `min` fd (p-1, q-1))
+    -- we either move on the q side, on the p side, or both
+
+--------------------------------------------------------------------------------
+-- * Dynamic programming / Memoization stuff
+
+-- | Create a table
+tabulate       :: Ix i => (i,i) -> (i -> a) -> Array i a
+tabulate rng f = listArray rng (map f $ range rng)
+
+-- | Memoize a function using an Array
+memo     :: Ix i => (i,i) -> (i -> a) -> (i -> a)
+memo rng = (!) . tabulate rng
+
+-- see: for the idea of this memoization
+-- https://byorgey.wordpress.com/2023/06/06/dynamic-programming-in-haskell-automatic-memoization/

hgeometry/test/PolyLine/Frechet/DiscreteSpec.hs

@@ -0,0 +1,29 @@
+module PolyLine.Frechet.DiscreteSpec where
+
+import qualified Data.List.NonEmpty as NonEmpty
+import           HGeometry.PolyLine
+import           HGeometry.PolyLine.Frechet.Discrete
+import           HGeometry.Point
+import           Test.Hspec
+
+--------------------------------------------------------------------------------
+
+type R = Int
+
+spec :: Spec
+spec = do
+         it "trivial" $ do
+           let ta = polyLineFromPoints' $ [origin, Point2 5 5]
+               tb = polyLineFromPoints'  $ [Point2 0 1, Point2 5 6, Point2 6 6]
+           discreteFrechetDistance ta tb `shouldBe` 2
+         it "trivial 2" $ do
+           let ta = polyLineFromPoints' [Point2 x 0     | x <- [1..10]]
+               tb = polyLineFromPoints' [Point2 x (y x) | x <- [1..10]]
+               y x = if x == 5 then 5 else 1
+           discreteFrechetDistance ta tb `shouldBe` 25
+
+discreteFrechetDistance :: PolyLine (Point 2 R) -> PolyLine (Point 2 R) -> R
+discreteFrechetDistance = frechetDistanceWith squaredEuclideanDist
+
+polyLineFromPoints' :: [Point 2 R] -> PolyLine (Point 2 R)
+polyLineFromPoints' = polyLineFromPoints . NonEmpty.fromList