Notebook source code: notebooks/10_practical_methods__shape_analysis.ipynb
Run it yourself on binder

Shape analysis of curves with the Square Root Velocity metric#

Lead author: Alice Le Brigant.

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 <geomstats/geomstats>`__ module. Translation and rotation can also be quotiented out using the align method of the `pre-shape.py <geomstats/geomstats>`__ 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.

 In [1]:

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

Example 1: plane curves#

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

 In [2]:

r2 = Euclidean(dim=2)
k_sampling_points = 20

curves_r2 = DiscreteCurves(ambient_manifold=r2, k_sampling_points=k_sampling_points)

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 [3]:

sampling_points = gs.linspace(0.0, 1.0, k_sampling_points)
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()

../_images/notebooks_10_practical_methods__shape_analysis_7_0.png

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 [4]:

curves_r2.metric.dist(point_a=curve_a, point_b=curve_b)

 Out [4]:

1.972134949925512

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

 In [5]:

geod_fun = curves_r2.metric.geodesic(curve_a, curve_b)

n_times = 20
times = gs.linspace(0.0, 1.0, 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()

../_images/notebooks_10_practical_methods__shape_analysis_12_0.png

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 [6]:

curve_a_resampled = parametrized_curve_a(sampling_points**2)
curve_b_resampled = parametrized_curve_b(sampling_points**2)

curves_r2.metric.dist(curve_a_resampled, curve_b_resampled)

 Out [6]:

1.9694017211949013

The geodesic keeps the same shape.

 In [7]:

geod_fun_1 = curves_r2.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()

../_images/notebooks_10_practical_methods__shape_analysis_16_0.png

However, if the curves are reparametrized in different ways, the distance changes, and so does the shape of the geodesic.

 In [8]:

curves_r2.metric.dist(curve_a, curve_b_resampled)

 Out [8]:

2.018033236326409

 In [9]:

geod_fun_2 = curves_r2.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()

../_images/notebooks_10_practical_methods__shape_analysis_19_0.png

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 [10]:

geod_velocity = n_times * (geod[1:] - geod[:-1])

curves_r2.equip_with_group_action("reparametrizations")
curves_r2.equip_with_quotient_structure()

geod_velocity_hor = curves_r2.fiber_bundle.horizontal_projection(geod_velocity, geod[:-1])
geod_velocity_ver = curves_r2.fiber_bundle.vertical_projection(geod_velocity, geod[:-1])

shooted_geod_fun = curves_r2.metric.geodesic(
    initial_point=curve_a, initial_tangent_vec=geod_velocity_hor[0]
)
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 [11]:

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

../_images/notebooks_10_practical_methods__shape_analysis_25_0.png

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 [12]:

hgeod_fun = curves_r2.metric.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()

../_images/notebooks_10_practical_methods__shape_analysis_27_0.png

We can check the horizontality of this geodesic by computing the norm of the vertical part of its velocity for all times.

 In [13]:

geod_vertical_norm = curves_r2.metric.norm(
    geod_velocity_ver, geod[:-1]
)

hgeod_velocity = n_times * (hgeod[1:] - hgeod[:-1])
hgeod_velocity_ver = curves_r2.fiber_bundle.vertical_projection(hgeod_velocity, hgeod[:-1])
hgeod_vertical_norm = curves_r2.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()

../_images/notebooks_10_practical_methods__shape_analysis_29_0.png

We can also check that this horizontal geodesic does not change if we resample the end curve.

 In [14]:

hgeod_fun = curves_r2.quotient.metric.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()

../_images/notebooks_10_practical_methods__shape_analysis_31_0.png

Finally, we can check that the quotient distance remains approximately constant for any parametrizations of the curves.

 In [15]:

print(curves_r2.quotient.metric.dist(curve_a, curve_b))
print(curves_r2.quotient.metric.dist(curve_a_resampled, curve_b))
print(curves_r2.quotient.metric.dist(curve_a_resampled, curve_b_resampled))
print(curves_r2.quotient.metric.dist(curve_a, curve_b_resampled))

1.7187665990832923
1.7195167839345487
1.7187052810816863
1.7178456824338564

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 [16]:

k_sampling_points = 100

r3 = Euclidean(dim=3)
curves_r3 = DiscreteCurves(ambient_manifold=r3, k_sampling_points=k_sampling_points)

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))
)

sampling_points = gs.linspace(0.0, 1.0, k_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()

../_images/notebooks_10_practical_methods__shape_analysis_36_0.png

 In [17]:

geod_fun = curves_r3.metric.geodesic(curve_a, curve_b)

n_times = 20
t = gs.linspace(0.0, 1.0, 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(k_sampling_points):
    ax.plot(geod[:, j, 0], geod[:, j, 1], geod[:, j, 2], "--", c="k")
plt.title("SRV geodesic")
plt.show()

../_images/notebooks_10_practical_methods__shape_analysis_37_0.png

 In [18]:

curves_r3.equip_with_group_action("reparametrizations")
curves_r3.equip_with_quotient_structure()

hgeod_fun = curves_r3.quotient.metric.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(k_sampling_points):
    ax.plot(hgeod[:, j, 0], hgeod[:, j, 1], hgeod[:, j, 2], "--", c="k")
plt.title("Horizontal SRV geodesic")
plt.show()

../_images/notebooks_10_practical_methods__shape_analysis_38_0.png

 In [19]:

geod_velocity = n_times * (geod[1:] - geod[:-1])
geod_velocity_ver = curves_r3.fiber_bundle.vertical_projection(geod_velocity, geod[:-1])
geod_vertical_norm = curves_r2.metric.norm(
    geod_velocity_ver, geod[:-1]
)

hgeod_velocity = n_times * (hgeod[1:] - hgeod[:-1])
hgeod_velocity_ver = curves_r3.fiber_bundle.vertical_projection(hgeod_velocity, hgeod[:-1])
hgeod_vertical_norm = curves_r2.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()

../_images/notebooks_10_practical_methods__shape_analysis_39_0.png

References#

[SKJJ2011]

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.

[L2017]

Le Brigant, “A discrete framework to find the optimal matching between manifold-valued curves,” in Journal of Mathematical Imaging and Vision, 61, pp. 40-70, 2019.