Tutorial: Shape analysis of curves with the Square Root Velocity metric

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 <https://github.com/geomstats/geomstats/blob/master/geomstats/geometry/discrete_curves.py>__ module. Translation and rotation can also be quotiented out using the align method of the pre-shape.py <https://github.com/geomstats/geomstats/blob/master/geomstats/geometry/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.

Setup

In :
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())
In :
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

Example 1: plane curves

We start with a basic example in $$\mathbb R^2$$.

In :
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.

In :
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() Distance between parametrized curves

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.

In :
curves_r2.square_root_velocity_metric.dist(point_a=curve_a, point_b=curve_b)
Out :
1.924866967405097

The distance, as any riemannian distance, is computed as the length of the geodesic.

In :
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.

In :
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)
Out :
1.9224597300658097

The geodesic keeps the same shape.

In :
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.

In :
curves_r2.square_root_velocity_metric.dist(curve_a, curve_b_resampled)
Out :
1.9696189997648135
In :
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() Distance between unparametrized curves

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 infinitesimal deformation of a curve can be split into the sum of vertical deformation (tangent to the fiber) that simply reparametrizes the curve without changing its shape, and a horizontal deformation (orthogonal to the fiber) that changes the shape. The distance between two unparametrized curves is then computed as the length of a horizontal geodesic linking their two fibers.

In practice, to compute the horizontal geodesic linking the fibers of two discrete parametrized curves curve_a and curve_b, we can fix the parametrization of curve_a, and search for a reparametrization of curve_b (i.e. another discrete curve with same shape as curve_b) that best “matches” curve_a.

Since geodesics that start with a horizontal velocity stay horizontal, a first idea would be the following:

• compute the geodesic between curve_a and curve_b

• compute the horizontal part of its initial velocity vector

• shoot from curve_a using this horizontal vector.

In :
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])
)

shooted_geod_fun = curves_r2.square_root_velocity_metric.geodesic(
initial_tangent_vec=geod_velocity_hor,
initial_curve=curve_a
)
shooted_geod = shooted_geod_fun(times)

The problem with this idea is that, while it yields a horizontal geodesic starting at curve_a, its end point does not belong to the fiber of curve_b: as we cas see below, the end curve of the horizontal geodesic is not a reparametrization of the initial red curve, it does not have the same shape.

In :
plt.figure(figsize=(7, 7))
plt.plot(shooted_geod[0, :, 0], shooted_geod[0, :, 1], 'o-b')
plt.plot(shooted_geod[-1, :, 0], shooted_geod[-1, :, 1], 'o-r')
for i in range(1, n_times - 1):
plt.plot(shooted_geod[i, :, 0], shooted_geod[i, :, 1], 'o-', c='k')
plt.title('Shooting with horizontal part of initial geodesic velocity vector')
plt.show() To obtain a horizontal geodesic starting at curve_a and ending at a curve with same shape as curve_b, we use an iterative optimal matching algorithm LAB2017. This algorithm moves along the fiber of curve_b to find the best representative with respect to curve_a by iterating the following steps:

• step 1: compute the geodesic between curve_a and the current representative of curve_b (initially, curve_b)

• step 2: compute the path whose velocity is a reparametrization of the horizontal part of the geodesic velocity at all time, and set the new representative of curve_b to be the end point of this path.

Note that the first step yields a geodesic that is not horizontal, while the second step yields a horizontal path that is not geodesic. By iterating these two steps, the algorithm converges to a horizontal geodesic.

In :
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.

In :
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.

In :
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.

In :
print(curves_r2.quotient_square_root_velocity_metric.dist(curve_a, curve_b))
print(curves_r2.quotient_square_root_velocity_metric.dist(curve_a_resampled, curve_b))
print(curves_r2.quotient_square_root_velocity_metric.dist(curve_a_resampled, curve_b_resampled))
print(curves_r2.quotient_square_root_velocity_metric.dist(curve_a, curve_b_resampled))
1.677393759010642
1.678054913356192
1.6774170579519256
1.6765232522160183

Example 2: 3D curves

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.

In :
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.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() In :
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.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() In :
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.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() In :
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() 