In this notebook, we demonstrate how to compute distances between curves in a way that does not depend on parametrization, i.e. that only depends on the shapes of the curves. This is achieved using the Square Root Velocity metric (see SKJJ2011) on the space of parametrized curves, and by quotienting out the action of reparametrization through an optimal matching algorithm (see LAB2017). We will use the discrete_curves.py
module. Translation and rotation can also be quotiented out using the align
method of the pre-shape.py
module, but we will not deal with these aspects here. See this usecase for details on the pre_shape.py
module, or this other usecase for an application where both modules are used.
import os
import subprocess
geomstats_gitroot_path = subprocess.check_output(
['git', 'rev-parse', '--show-toplevel'],
universal_newlines=True)
os.chdir(geomstats_gitroot_path[:-1])
print('Working directory: ', os.getcwd())
Working directory: /Users/alicelebrigant/ownCloud/Python/SpyderProjects/geomstats
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import geomstats.backend as gs
from geomstats.geometry.euclidean import Euclidean
from geomstats.geometry.discrete_curves import DiscreteCurves
INFO: Using numpy backend
We start with a basic example in $\mathbb R^2$.
r2 = Euclidean(dim=2)
curves_r2 = DiscreteCurves(ambient_manifold=r2)
parametrized_curve_a = lambda x: gs.transpose(gs.array([1 + 2 * gs.sin(gs.pi * x), 3 + 2 * gs.cos(gs.pi * x)]))
parametrized_curve_b = lambda x: gs.transpose(gs.array([5 * gs.ones(len(x)), 4 * (1 - x) + 1]))
In practice, we work with discrete curves, i.e. sample points from the parametrized curves.
n_sampling_points = 20
sampling_points = gs.linspace(0., 1., n_sampling_points + 1)
curve_a = parametrized_curve_a(sampling_points)
curve_b = parametrized_curve_b(sampling_points)
plt.figure(figsize=(7, 7))
plt.plot(curve_a[:, 0], curve_a[:, 1], 'o-b')
plt.plot(curve_b[:, 0], curve_b[:, 1], 'o-r')
plt.show()
The metric we use to compare parametrized curves is the so-called Square Root Velocity metric, that computes an $L^2$ distance between the velocities of the curves, suitably renormalized to get reparametrization invariance. See SKJJ2011 for more details.
curves_r2.square_root_velocity_metric.dist(point_a=curve_a, point_b=curve_b)
array([9.68542047])
The distance, as any riemannian distance, is computed as the length of the geodesic.
geod_fun = curves_r2.square_root_velocity_metric.geodesic(initial_curve=curve_a, end_curve=curve_b)
n_times = 20
times = gs.linspace(0., 1., n_times)
geod = geod_fun(times)
plt.figure(figsize=(7, 7))
plt.plot(geod[0, :, 0], geod[0, :, 1], 'o-b')
for i in range(1, n_times - 1):
plt.plot(geod[i, :, 0], geod[i, :, 1], 'o-k')
plt.plot(geod[-1, :, 0], geod[-1, :, 1], 'o-r')
plt.title('Geodesic for the Square Root Velocity metric')
plt.show()
The Square Root Velocity metric is reparametrization invariant in the sense that, if the two curves are reparametrized in the same way, the distance does not change.
curve_a_resampled = parametrized_curve_a(sampling_points ** 2)
curve_b_resampled = parametrized_curve_b(sampling_points ** 2)
curves_r2.square_root_velocity_metric.dist(curve_a_resampled, curve_b_resampled)
array([9.67537491])
The geodesic keeps the same shape.
geod_fun_1 = curves_r2.square_root_velocity_metric.geodesic(curve_a_resampled, curve_b_resampled)
geod_1 = geod_fun_1(times)
plt.figure(figsize=(7, 7))
plt.plot(geod_1[0, :, 0], geod_1[0, :, 1], 'o-b')
for i in range(1, n_times - 1):
plt.plot(geod_1[i, :, 0], geod_1[i, :, 1], 'o-k')
plt.plot(geod_1[-1, :, 0], geod_1[-1, :, 1], 'o-r')
plt.title('Geodesic when both curves are reparametrized in the same way')
plt.show()
However, if the curves are reparametrized in different ways, the distance changes, and so does the shape of the geodesic.
curves_r2.square_root_velocity_metric.dist(curve_a, curve_b_resampled)
array([9.87255687])
geod_fun_2 = curves_r2.square_root_velocity_metric.geodesic(curve_a, curve_b_resampled)
geod_2 = geod_fun_2(times)
plt.figure(figsize=(7, 7))
plt.plot(geod_2[0, :, 0], geod_2[0, :, 1], 'o-b')
for i in range(1, n_times - 1):
plt.plot(geod_2[i, :, 0], geod_2[i, :, 1], 'o-k')
plt.plot(geod_2[-1, :, 0], geod_2[-1, :, 1], 'o-r')
plt.title('Geodesic when only the red curve is reparametrized')
plt.show()
In order to completely quotient out parametrization, distances are computed in the base space of a fiber bundle where the fibers represent equivalent classes of curves with the same shape (i.e. equal modulo reparametrization). Any tangent vector to a curve (i.e. infinitesimal vector field along that curve) can be split into the sum of vertical vector, tangent to the fiber, and a horizontal vector, orthogonal to all vertical vectors. The distance between two unparametrized curves is then computed as the length of a horizontal geodesic linking their two fibers. In practice, given two discrete versions of parametrized curves curve_a
and curve_b
, we search for the horizontal geodesic that starts from curve_a
and that ends at a curve that belongs to the same fiber as curve_b
. This horizontal geodesic is computed iteratively by an optimal matching algorithm, see LAB2017.
hgeod_fun = curves_r2.quotient_square_root_velocity_metric.horizontal_geodesic(curve_a, curve_b)
hgeod = hgeod_fun(times)
plt.figure(figsize=(7, 7))
plt.plot(hgeod[0, :, 0], hgeod[0, :, 1], 'o-b')
for i in range(1, n_times - 1):
plt.plot(hgeod[i, :, 0], hgeod[i, :, 1], 'o-k')
plt.plot(hgeod[-1, :, 0], hgeod[-1, :, 1], 'o-r')
plt.title('Horizontal geodesic')
plt.show()
We can check the horizontality of this geodesic by computing the norm of the vertical part of its velocity for all times.
geod_velocity = n_times * (geod[1:] - geod[:-1])
geod_velocity_hor, geod_velocity_ver, _ = (
curves_r2.quotient_square_root_velocity_metric.split_horizontal_vertical(geod_velocity, geod[:-1])
)
geod_vertical_norm = curves_r2.square_root_velocity_metric.norm(geod_velocity_ver, geod[:-1])
hgeod_velocity = n_times * (hgeod[1:] - hgeod[:-1])
hgeod_velocity_hor, hgeod_velocity_ver, _ = (
curves_r2.quotient_square_root_velocity_metric.split_horizontal_vertical(hgeod_velocity, hgeod[:-1])
)
hgeod_vertical_norm = curves_r2.square_root_velocity_metric.norm(hgeod_velocity_ver, hgeod[:-1])
plt.figure()
plt.plot(times[:-1], geod_vertical_norm, 'o', label='initial geodesic')
plt.plot(times[:-1], hgeod_vertical_norm, 'o', label='horizontal geodesic')
plt.legend()
plt.title('Norm of the vertical part of the geodesic velocity')
plt.show()
We can also check that this horizontal geodesic does not change if we resample the end curve.
hgeod_fun = curves_r2.quotient_square_root_velocity_metric.horizontal_geodesic(curve_a, curve_b_resampled)
hgeod = hgeod_fun(times)
plt.figure(figsize=(7, 7))
plt.plot(hgeod[0, :, 0], hgeod[0, :, 1], 'o-b')
for i in range(1, n_times - 1):
plt.plot(hgeod[i, :, 0], hgeod[i, :, 1], 'o-k')
plt.plot(hgeod[-1, :, 0], hgeod[-1, :, 1], 'o-r')
plt.title('Horizontal geodesic when the red curve is reparametrized')
plt.show()
Finally, we can check that the quotient distance remains approximately constant for any parametrizations of the curves.
print(curves_r2.quotient_square_root_velocity_metric.quotient_dist(curve_a, curve_b))
print(curves_r2.quotient_square_root_velocity_metric.quotient_dist(curve_a_resampled, curve_b))
print(curves_r2.quotient_square_root_velocity_metric.quotient_dist(curve_a_resampled, curve_b_resampled))
print(curves_r2.quotient_square_root_velocity_metric.quotient_dist(curve_a, curve_b_resampled))
1.718817126370995 1.71949460797375 1.718841000681418 1.717925122344401
Below we follow similar steps for curves in $\mathbb R^3$. In this example, we can see that the horizontal geodesic "straightens out" the original geodesic.
r3 = Euclidean(dim=3)
curves_r3 = DiscreteCurves(ambient_manifold=r3)
parametrized_curve_a = lambda x: gs.transpose(gs.stack((gs.cos(2 + 8 * x), gs.sin(2 + 8 * x), 2 + 10 * x)))
parametrized_curve_b = lambda x: gs.transpose(gs.stack((gs.cos(4 + 8 * x), gs.sin(4 + 8 * x), 2 + 10 * x)))
n_sampling_points = 100
sampling_points = gs.linspace(0., 1., n_sampling_points)
curve_a = parametrized_curve_a(sampling_points)
curve_b = parametrized_curve_b(sampling_points)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.plot(curve_a[:, 0], curve_a[:, 1], curve_a[:, 2], 'b')
ax.plot(curve_b[:, 0], curve_b[:, 1], curve_b[:, 2], 'r')
ax.scatter(curve_a[0, 0], curve_a[0, 1], curve_a[0, 2], 'b')
ax.scatter(curve_b[0, 0], curve_b[0, 1], curve_b[0, 2], 'r')
plt.show()
geod_fun = curves_r3.square_root_velocity_metric.geodesic(initial_curve=curve_a, end_curve=curve_b)
n_times = 20
t = gs.linspace(0., 1., n_times)
geod = geod_fun(t)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.plot(curve_a[:, 0], curve_a[:, 1], curve_a[:, 2], '-', c='b', linewidth=2)
ax.plot(curve_b[:, 0], curve_b[:, 1], curve_b[:, 2], '-', c='r', linewidth=2)
for i in range(1, n_times - 1):
ax.plot(geod[i, :, 0], geod[i, :, 1], geod[i, :, 2], '-', c='k')
for j in range(n_sampling_points):
ax.plot(geod[:, j, 0], geod[:, j, 1], geod[:, j, 2], '--', c='k')
plt.title('SRV geodesic')
plt.show()
hgeod_fun = curves_r3.quotient_square_root_velocity_metric.horizontal_geodesic(curve_a, curve_b)
hgeod = hgeod_fun(t)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.plot(curve_a[:, 0], curve_a[:, 1], curve_a[:, 2], '-', c='b', linewidth=2)
ax.plot(curve_b[:, 0], curve_b[:, 1], curve_b[:, 2], '-', c='r', linewidth=2)
for i in range(1, n_times - 1):
ax.plot(hgeod[i, :, 0], hgeod[i, :, 1], hgeod[i, :, 2], '-', c='k')
for j in range(n_sampling_points):
ax.plot(hgeod[:, j, 0], hgeod[:, j, 1], hgeod[:, j, 2], '--', c='k')
plt.title('Horizontal SRV geodesic')
plt.show()
geod_velocity = n_times * (geod[1:] - geod[:-1])
geod_velocity_hor, geod_velocity_ver, _ = (
curves_r2.quotient_square_root_velocity_metric.split_horizontal_vertical(geod_velocity, geod[:-1])
)
geod_vertical_norm = curves_r2.square_root_velocity_metric.norm(geod_velocity_ver, geod[:-1])
hgeod_velocity = n_times * (hgeod[1:] - hgeod[:-1])
hgeod_velocity_hor, hgeod_velocity_ver, _ = (
curves_r2.quotient_square_root_velocity_metric.split_horizontal_vertical(hgeod_velocity, hgeod[:-1])
)
hgeod_vertical_norm = curves_r2.square_root_velocity_metric.norm(hgeod_velocity_ver, hgeod[:-1])
plt.figure()
plt.plot(times[:-1], geod_vertical_norm, 'o', label='initial geodesic')
plt.plot(times[:-1], hgeod_vertical_norm, 'o', label='horizontal geodesic')
plt.legend()
plt.title('Norm of the vertical part of the geodesic velocity')
plt.show()
.. [SKJJ2011] A. Srivastava, E. Klassen, S. H. Joshi and I. H. Jermyn, "Shape Analysis of Elastic Curves in Euclidean Spaces," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 33, no. 7, pp. 1415-1428, July 2011.
.. [LAB2017] A. Le Brigant, M. Arnaudon and F. Barbaresco, "Optimal matching between curves in a manifold," in International Conference on Geometric Science of Information, pp. 57-65, Springer, Cham, 2017.