弗雷歇距离的原理及python代码实现（动态规划）

弗雷歇距离的原理及python代码实现（动态规划）

在网上看了很多关于弗雷歇距离的介绍，结合自己的理解，出一版更通俗易懂、更清晰具体的解释。

最简单的解释自然是最短狗绳长度，但我将从另一个角度来解释它。

图中人牵着狗在走，人走直线，狗走得自由且散漫。为了能拴住狗，任何时刻狗绳的长度都应该大于人狗距离，于是有最短狗绳长度等于最大人狗距离。

现在我们假定人只能走蓝色轨道，狗只能走红色轨道且都只能向前走，但是具体怎么走——中途停不停、走多快是未知的，然后我们对时间进行采样，将得到人、狗轨迹的离散序列，计算两序列中对应离散点的距离，距离的最大值即这一次遛狗所需的最短狗绳长度。

基于这个理解，上图显示的便是一次遛狗所得到的人、狗轨迹离散序列，虚线连接了两序列中对应的离散点，最大虚线长度便是这次遛狗所需的最短狗绳长度。下一次遛狗得到的离散序列可能就大不相同，例如，人走到中间时，狗已经走到红色轨道的终点并在那等着主人，那么这一次遛狗所需的最短狗绳长度就大于上一次的。在无数次遛狗中，最短狗绳长度的最小值就是弗雷歇距离！

注：对于两条连续轨迹，不可能遍历所有人狗轨迹离散序列，因此只能求下界，对于两条离散轨迹，取最小值即可。

动态规划计算弗雷歇距离：

设甲、乙两人的离散轨迹序列分别为Ta、Tb，Ta、Tb的长度分别为M、N，两轨迹之间的弗雷歇距离为Fd[M][N]。假设Ta、Tb的最后一个点分别是A、B，那么可以有三种行走情形：

1. 甲、乙两人都在最后时刻到达终点A、B，那么此时Fd[M][N] = max(Fd[M-1][N-1], d(A,B))
2. 甲先到达终点A并在那等着乙到达终点B，那么此时Fd[M][N] = max(Fd[M][N-1], d(A,B))
3. 乙先到达终点B并在那等着甲到达终点A，那么此时Fd[M][N] = max(Fd[M-1][N], d(A,B))

由此我们可以得到传递函数：
$$Fd[i][j]=\{\begin{array}{lcl}d(Ta[i],Tb[j])& & i=j=1\\ max(Fd[i][j-1],d(Ta[i],Tb[j]))& & i=1,j>1\\ max(Fd[i-1][j],d(Ta[i],Tb[j]))& & i>1,j=1\\ max(min(Fd[i-1][j-1],Fd[i-1][j],Fd[i][j-1]),d(Ta[i],Tb[j]))& & i>1,j>1\end{array}$$

$$\begin{array}{l}min(max(Fd[i-1][j-1],d(Ta[i],Tb[j])),max(Fd[i-1][j],d(Ta[i],Tb[j])),max(Fd[i][j-1],d(Ta[i],Tb[j])))\\ =max(min(Fd[i-1][j-1],Fd[i-1][j],Fd[i][j-1]),d(Ta[i],Tb[j]))\end{array}$$

python代码实现(https://blog.csdn.net/qq_42517334/article/details/103506868)：

import math
import numpy as np


def calculate_euclid(point_a, point_b):
    """
    Args:
        point_a: a data point of curve_a
        point_b: a data point of curve_b
    Return:
        The Euclid distance between point_a and point_b
    """
    return math.sqrt((point_a - point_b) ** 2)


def calculate_frechet_distance(dp, i, j, curve_a, curve_b):
    """
    Args:
        dp: The distance matrix
        i: The index of curve_a
        j: The index of curve_b
        curve_a: The data sequence of curve_a
        curve_b: The data sequence of curve_b
    Return:
        The frechet distance between curve_a[i] and curve_b[j]
    """
    if dp[i][j] > -1:
        return dp[i][j]
    elif i == 0 and j == 0:
        dp[i][j] = calculate_euclid(curve_a[0], curve_b[0])
    elif i > 0 and j == 0:
        dp[i][j] = max(calculate_frechet_distance(dp, i - 1, 0, curve_a, curve_b),
                       calculate_euclid(curve_a[i], curve_b[0]))
    elif i == 0 and j > 0:
        dp[i][j] = max(calculate_frechet_distance(dp, 0, j - 1, curve_a, curve_b),
                       calculate_euclid(curve_a[0], curve_b[j]))
    elif i > 0 and j > 0:
        dp[i][j] = max(min(calculate_frechet_distance(dp, i - 1, j, curve_a, curve_b),
                           calculate_frechet_distance(dp, i - 1, j - 1, curve_a, curve_b),
                           calculate_frechet_distance(dp, i, j - 1, curve_a, curve_b)),
                       calculate_euclid(curve_a[i], curve_b[j]))
    else:
        dp[i][j] = float("inf")
    return dp[i][j]


def get_similarity(curve_a, curve_b):
    dp = [[-1 for _ in range(len(curve_b))] for _ in range(len(curve_a))]
    similarity = calculate_frechet_distance(dp, len(curve_a) - 1, len(curve_b) - 1, curve_a, curve_b)
    # return max(np.array(dp).reshape(-1, 1))[0]
    return similarity


if __name__ == '__main__':
    Ta = [1, 2, 3, 4]
    Tb = [7, 8, 9]
    print(get_similarity(Ta, Tb))

Reference：

数学定义

代码来源