In [1]:
# !/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on 20181221
@author: zhangji
"""
%pylab inline
# pylab.rcParams['figure.figsize'] = (25, 11)
fontsize = 40
import numpy as np
import math
import scipy as sp
from scipy.optimize import leastsq, curve_fit
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy.io import loadmat, savemat
# import scipy.misc
import matplotlib
from matplotlib import pyplot as plt
from matplotlib import animation, rc
import matplotlib.ticker as mtick
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
from mpl_toolkits.mplot3d import Axes3D, axes3d
from sympy import symbols, simplify, series, exp
from sympy.matrices import Matrix
from sympy.solvers import solve
from IPython.display import display, HTML
from tqdm import tqdm_notebook as tqdm
import pandas as pd
import re
from scanf import scanf
import os
import glob
import natsort
from shutil import copyfile
from codeStore import support_fun as spf
from src.support_class import *
from src import stokes_flow as sf
rc('animation', html='html5')
PWD = os.getcwd()
font = {'size': 20}
matplotlib.rc('font', **font)
np.set_printoptions(linewidth=90, precision=5)
In [63]:
def _c(ca,i,j,p,q):
if ca[i,j] > -1:
return ca[i,j]
elif i == 0 and j == 0:
ca[i,j] = np.linalg.norm(p[i]-q[j])
elif i > 0 and j == 0:
ca[i,j] = max( _c(ca,i-1,0,p,q), np.linalg.norm(p[i]-q[j]) )
elif i == 0 and j > 0:
ca[i,j] = max( _c(ca,0,j-1,p,q), np.linalg.norm(p[i]-q[j]) )
elif i > 0 and j > 0:
ca[i,j] = max( \
min( \
_c(ca,i-1,j,p,q), \
_c(ca,i-1,j-1,p,q), \
_c(ca,i,j-1,p,q) \
), \
np.linalg.norm(p[i]-q[j]) \
)
else:
ca[i,j] = float('inf')
return ca[i,j]
def frdist(p,q):
"""
Computes the discrete Fréchet distance between
two curves. The Fréchet distance between two curves in a
metric space is a measure of the similarity between the curves.
The discrete Fréchet distance may be used for approximately computing
the Fréchet distance between two arbitrary curves,
as an alternative to using the exact Fréchet distance between a polygonal
approximation of the curves or an approximation of this value.
This is a Python 3.* implementation of the algorithm produced
in Eiter, T. and Mannila, H., 1994. Computing discrete Fréchet distance. Tech.
Report CD-TR 94/64, Information Systems Department, Technical University
of Vienna.
http://www.kr.tuwien.ac.at/staff/eiter/et-archive/cdtr9464.pdf
Function dF(P, Q): real;
input: polygonal curves P = (u1, . . . , up) and Q = (v1, . . . , vq).
return: δdF (P, Q)
ca : array [1..p, 1..q] of real;
function c(i, j): real;
begin
if ca(i, j) > −1 then return ca(i, j)
elsif i = 1 and j = 1 then ca(i, j) := d(u1, v1)
elsif i > 1 and j = 1 then ca(i, j) := max{ c(i − 1, 1), d(ui, v1) }
elsif i = 1 and j > 1 then ca(i, j) := max{ c(1, j − 1), d(u1, vj ) }
elsif i > 1 and j > 1 then ca(i, j) :=
max{ min(c(i − 1, j), c(i − 1, j − 1), c(i, j − 1)), d(ui, vj ) }
else ca(i, j) = ∞
return ca(i, j);
end; /* function c */
begin
for i = 1 to p do for j = 1 to q do ca(i, j) := −1.0;
return c(p, q);
end.
Parameters
----------
P : Input curve - two dimensional array of points
Q : Input curve - two dimensional array of points
Returns
-------
dist: float64
The discrete Fréchet distance between curves `P` and `Q`.
Examples
--------
>>> from frechetdist import frdist
>>> P=[[1,1], [2,1], [2,2]]
>>> Q=[[2,2], [0,1], [2,4]]
>>> frdist(P,Q)
>>> 2.0
>>> P=[[1,1], [2,1], [2,2]]
>>> Q=[[1,1], [2,1], [2,2]]
>>> frdist(P,Q)
>>> 0
"""
p = np.array(p, np.float64)
q = np.array(q, np.float64)
len_p = len(p)
len_q = len(q)
if len_p == 0 or len_q == 0:
raise ValueError('Input curves are empty.')
if len_p != len_q or len(p[0]) != len(q[0]):
raise ValueError('Input curves do not have the same dimensions.')
ca = ( np.ones((len_p,len_q), dtype=np.float64) * -1 )
dist = _c(ca,len_p-1,len_q-1,p,q)
return dist
# Euclidean distance.
def euc_dist(pt1,pt2):
return math.sqrt((pt2[0]-pt1[0])*(pt2[0]-pt1[0])+(pt2[1]-pt1[1])*(pt2[1]-pt1[1]))
def _c(ca,i,j,P,Q):
if ca[i,j] > -1:
return ca[i,j]
elif i == 0 and j == 0:
ca[i,j] = euc_dist(P[0],Q[0])
elif i > 0 and j == 0:
ca[i,j] = max(_c(ca,i-1,0,P,Q),euc_dist(P[i],Q[0]))
elif i == 0 and j > 0:
ca[i,j] = max(_c(ca,0,j-1,P,Q),euc_dist(P[0],Q[j]))
elif i > 0 and j > 0:
ca[i,j] = max(min(_c(ca,i-1,j,P,Q),_c(ca,i-1,j-1,P,Q),_c(ca,i,j-1,P,Q)),euc_dist(P[i],Q[j]))
else:
ca[i,j] = float("inf")
return ca[i,j]
""" Computes the discrete frechet distance between two polygonal lines
Algorithm: http://www.kr.tuwien.ac.at/staff/eiter/et-archive/cdtr9464.pdf
P and Q are arrays of 2-element arrays (points)
"""
def frechetDist(P,Q):
ca = np.ones((len(P),len(Q)))
ca = np.multiply(ca,-1)
return _c(ca,len(P)-1,len(Q)-1,P,Q)
def read_ecoli_mat(mat_name):
mat_contents = loadmat(mat_name)
ecoli_U = mat_contents['ecoli_U']
ecoli_norm = mat_contents['ecoli_norm']
ecoli_center = mat_contents['ecoli_center']
return ecoli_center, ecoli_norm, ecoli_U
In [59]:
P=[[1,1], [2,1], [2,2]]
Q=[[1,1], [2,1], [2,2], [2,2]]
frechetDist(P,Q)
Out[59]:
In [73]:
base_mat = os.path.join(PWD, 'ecoli_shear1c', 'eq_dt0.010_O5.mat')
dir_name = 'ecoli_shear1c'
base_center, base_norm, base_U = read_ecoli_mat(base_mat)
base_length = np.linalg.norm((base_center[:-1, :] - base_center[1:, :]), axis=1).sum()
_, dt0, _ = scanf('%s/eq_dt%f_%s', base_mat)
t_dir = os.path.join(PWD, dir_name)
mat_names = glob.glob('%s/*.mat' % t_dir)
for mati in natsort.natsorted(mat_names):
ecoli_center, ecoli_norm, ecoli_U = read_ecoli_mat(mati)
_, dt, _ = scanf('%s/eq_dt%f_%s', mati)
scale_cut = int(ecoli_center.shape[0] // (dt / dt0))
t_dst = frechetDist(ecoli_center[:scale_cut, :], base_center)
print(mati, t_dst, t_dst / base_length)
In [71]:
np.linalg.norm((base_center[:-1, :] - base_center[1:, :]), axis=1).sum()
Out[71]: