0 people like it.

Discrete Fréchet Distance

Compute the Discrete Fréchet Distance between two arrays (which may be of different lengths). Based on the 1994 algorithm by Thomas Eiter and Heikki Mannila. Not extensively tested, so use at your peril! (This version with some small fixes.)

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54: 
55: 
56: 
57: 
58: 
59: 
60: 
61: 
62: 
63: 
64: 
65: 

open System

// Discrete Frechet Distance: 
// (Based on the 1994 algorithm from Thomas Eiter and Heikki Mannila.)
let frechet (P : array<float*float>) (Q : array<float*float>) =
    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
namespace System
val frechet : P:(float * float) array -> Q:(float * float) array -> float

Full name: Script.frechet
val P : (float * float) array
type 'T array = 'T []

Full name: Microsoft.FSharp.Core.array<_>
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
val Q : (float * float) array
val sq : (float -> float)
val x : float
val min3 : ('a -> 'a -> 'a -> 'a) (requires comparison)
val x : 'a (requires comparison)
val y : 'a (requires comparison)
val z : 'a (requires comparison)
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
  interface IEnumerable
  interface IEnumerable<'T>
  member Head : 'T
  member IsEmpty : bool
  member Item : index:int -> 'T with get
  member Length : int
  member Tail : 'T list
  static member Cons : head:'T * tail:'T list -> 'T list
  static member Empty : 'T list

Full name: Microsoft.FSharp.Collections.List<_>
val min : list:'T list -> 'T (requires comparison)

Full name: Microsoft.FSharp.Collections.List.min
val d : (float * float -> float * float -> float)
val a : float * float
val b : float * float
val ab_x : float
val fst : tuple:('T1 * 'T2) -> 'T1

Full name: Microsoft.FSharp.Core.Operators.fst
val ab_y : float
val snd : tuple:('T1 * 'T2) -> 'T2

Full name: Microsoft.FSharp.Core.Operators.snd
val sqrt : value:'T -> 'U (requires member Sqrt)

Full name: Microsoft.FSharp.Core.Operators.sqrt
val p : int
val q : int
type Array =
  member Clone : unit -> obj
  member CopyTo : array:Array * index:int -> unit + 1 overload
  member GetEnumerator : unit -> IEnumerator
  member GetLength : dimension:int -> int
  member GetLongLength : dimension:int -> int64
  member GetLowerBound : dimension:int -> int
  member GetUpperBound : dimension:int -> int
  member GetValue : [<ParamArray>] indices:int[] -> obj + 7 overloads
  member Initialize : unit -> unit
  member IsFixedSize : bool
  ...

Full name: System.Array
val length : array:'T [] -> int

Full name: Microsoft.FSharp.Collections.Array.length
val ca : float [,]
module Array2D

from Microsoft.FSharp.Collections
val init : length1:int -> length2:int -> initializer:(int -> int -> 'T) -> 'T [,]

Full name: Microsoft.FSharp.Collections.Array2D.init
val c : (int -> int -> float)
val i : int
val j : int
type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math
Math.Max(val1: decimal, val2: decimal) : decimal
   (+0 other overloads)
Math.Max(val1: float, val2: float) : float
   (+0 other overloads)
Math.Max(val1: float32, val2: float32) : float32
   (+0 other overloads)
Math.Max(val1: uint64, val2: uint64) : uint64
   (+0 other overloads)
Math.Max(val1: int64, val2: int64) : int64
   (+0 other overloads)
Math.Max(val1: uint32, val2: uint32) : uint32
   (+0 other overloads)
Math.Max(val1: int, val2: int) : int
   (+0 other overloads)
Math.Max(val1: uint16, val2: uint16) : uint16
   (+0 other overloads)
Math.Max(val1: int16, val2: int16) : int16
   (+0 other overloads)
Math.Max(val1: byte, val2: byte) : byte
   (+0 other overloads)
val nan : float

Full name: Microsoft.FSharp.Core.Operators.nan
val a1 : (float * float) array
val a2 : (float * float) array
val abs : value:'T -> 'T (requires member Abs)

Full name: Microsoft.FSharp.Core.Operators.abs
val linearPositive : (float * float) []

Full name: Script.linearPositive
Multiple items
val seq : sequence:seq<'T> -> seq<'T>

Full name: Microsoft.FSharp.Core.Operators.seq

--------------------
type seq<'T> = Collections.Generic.IEnumerable<'T>

Full name: Microsoft.FSharp.Collections.seq<_>
val ofSeq : source:seq<'T> -> 'T []

Full name: Microsoft.FSharp.Collections.Array.ofSeq
val linearNegative : (float * float) []

Full name: Script.linearNegative
val linearPositiveOffset10 : (float * float) []

Full name: Script.linearPositiveOffset10
val x_x2 : (float * float) []

Full name: Script.x_x2
val x_x2_plus2y : (float * float) []

Full name: Script.x_x2_plus2y
val x_x2_plus2x : (float * float) []

Full name: Script.x_x2_plus2x
val x_x2_plus10xy : (float * float) []

Full name: Script.x_x2_plus10xy
val x_x2_scaled_x : (float * float) []

Full name: Script.x_x2_scaled_x
val circle : (float * float) []

Full name: Script.circle
field Math.PI = 3.14159265359
Math.Sin(a: float) : float
Math.Cos(d: float) : float
val circleStretched : (float * float) []

Full name: Script.circleStretched
val circleStretchedHiRes : (float * float) []

Full name: Script.circleStretchedHiRes
val test1 : float

Full name: Script.test1
val test2 : float

Full name: Script.test2
val test3 : float

Full name: Script.test3
val test4 : float

Full name: Script.test4
val test5 : float

Full name: Script.test5
val test6 : float

Full name: Script.test6
val test7 : float

Full name: Script.test7
val test8 : float

Full name: Script.test8
val test9 : float

Full name: Script.test9
val test10 : float

Full name: Script.test10

More information

Link:http://fssnip.net/am
Posted:12 years ago
Author:Kit Eason
Tags: frechet distance , curves , math

This web site is created using F# and Suave web server. It is hosted on Azure and the source code is on GitHub. Contributions are welcome!

The first version of fssnip.net has been created by @tomaspetricek back in 2010. This web site is a new, open-source and contribution-friendly version.