open System // Discrete Frechet Distance: // (Based on the 1994 algorithm from Thomas Eiter and Heikki Mannila.) let frechet (P : array) (Q : array) = let sq (x : float) = x * x let min3 x y z = [x; y; z] |> List.min let d (a : float*float) (b: float*float) = let ab_x = fst(a) - fst(b) let ab_y = snd(a) - snd(b) sqrt(sq ab_x + sq ab_y) let p, q = Array.length P, Array.length Q let ca = Array2D.init p q (fun _ _ -> -1.0) let rec c i j = if ca.[i, j] > -1.0 then ca.[i, j] else if i = 0 && j = 0 then ca.[i, j] <- d (P.[0]) (Q.[0]) elif i > 0 && j = 0 then ca.[i, j] <- Math.Max((c (i-1) 0), (d P.[i] Q.[0])) elif i = 0 && j > 0 then ca.[i, j] <- Math.Max((c 0 (j-1)), (d P.[0] Q.[j])) elif i > 0 && j > 0 then ca.[i, j] <- Math.Max(min3 (c (i-1) j) (c (i-1) (j-1)) (c i (j-1)), (d P.[i] Q.[j])) else ca.[i, j] <- nan ca.[i, j] c (p-1) (q-1) // Use frechet as an operator: let (-~~) a1 a2 = abs(frechet a1 a2) // Test arrays: let linearPositive = (seq {for x in 0. .. 1. .. 100. do yield (x, x)}) |> Array.ofSeq let linearNegative = (seq {for x in 0. .. 1. .. 100. do yield (x, 100.-x)}) |> Array.ofSeq let linearPositiveOffset10 = (seq {for x in 0. .. 1. .. 100. do yield (x+10., x+10.)}) |> Array.ofSeq let x_x2 = (seq {for x in 0. .. 1. .. 100. do yield (x, x*x)}) |> Array.ofSeq let x_x2_plus2y = (seq {for x in 0. .. 1. .. 100. do yield (x, x*x+2.)}) |> Array.ofSeq let x_x2_plus2x = (seq {for x in 0. .. 1. .. 100. do yield (x+2., x*x)}) |> Array.ofSeq let x_x2_plus10xy = (seq {for x in 0. .. 1. .. 100. do yield (x+10., x*x+10.)}) |> Array.ofSeq let x_x2_scaled_x = (seq {for x in 0. .. 1. .. 100. do yield (x/2., x*x+10.)}) |> Array.ofSeq let circle = (seq {for x in 0. .. 0.1 .. System.Math.PI*2. do yield (Math.Sin(x), Math.Cos(x))}) |> Array.ofSeq let circleStretched = (seq {for x in 0. .. 0.1 .. System.Math.PI*2. do yield (Math.Sin(x)*4.5, Math.Cos(x)*3.4)}) |> Array.ofSeq let circleStretchedHiRes = (seq {for x in 0. .. 0.01 .. System.Math.PI*2. do yield (Math.Sin(x)*4.5, Math.Cos(x)*3.4)}) |> Array.ofSeq // Tests: let test1 = linearPositive -~~ linearPositive // 0.0 - identical let test2 = linearPositive -~~ linearNegative // 100.0 let test3 = linearPositive -~~ linearPositiveOffset10 // 14.14213562 let test4 = x_x2 -~~ x_x2 // 0.0 - identical let test5 = x_x2 -~~ x_x2_plus2y // 2.0 let test6 = x_x2 -~~ x_x2_plus2x // 2.0 let test7 = x_x2 -~~ x_x2_plus10xy // 14.14213562 let test8 = x_x2 -~~ x_x2_scaled_x // 50.99019514 let test9 = circle -~~ circleStretched // 3.4998577 let test10 = circle -~~ circleStretchedHiRes // 3.500879039