Structured Deformations to Model Acropetal Growth

1. Curve and dots registration using implicit modules of order 1, learning the growth model tensor.

2. Curve registration using implicit modules of order with learned growth model tensor.

Import relevant Python modules.

import sys
sys.path.append("../../")
import math
import copy
import pickle

import torch

import matplotlib.pyplot as plt

import imodal

torch.set_default_dtype(torch.float32)
imodal.Utilities.set_compute_backend('torch')

Learning the growth model tensor

We load the data (shape and dots of the source and target leaves).

with open("../../data/acropetal.pickle", 'rb') as f:
    data = pickle.load(f)

dots_source = torch.tensor(data['dots_source'], dtype=torch.get_default_dtype())
dots_target = torch.tensor(data['dots_target'], dtype=torch.get_default_dtype())
shape_source = imodal.Utilities.close_shape(torch.tensor(data['shape_source']).type(torch.get_default_dtype()))
shape_target = imodal.Utilities.close_shape(torch.tensor(data['shape_target']).type(torch.get_default_dtype()))


aabb_source = imodal.Utilities.AABB.build_from_points(shape_source)
aabb_target = imodal.Utilities.AABB.build_from_points(shape_target)

Plot source and target.

plt.title("Source and target")
plt.plot(shape_source[:, 0].numpy(), shape_source[:, 1].numpy(), color='black')
plt.plot(dots_source[:, 0].numpy(), dots_source[:, 1].numpy(), '.', color='black')
plt.plot(shape_target[:, 0].numpy(), shape_target[:, 1].numpy(), color='red')
plt.plot(dots_target[:, 0].numpy(), dots_target[:, 1].numpy(), '.', color='red')

plt.axis('equal')
plt.show()

We now sample the points that will be used by the implicit deformation module of order 1 (growth module).

# Build AABB (Axis Aligned Bounding Box) around the source shape and uniformly
# sample points for the growth module.
points_density = 0.05

aabb_source = imodal.Utilities.AABB.build_from_points(shape_source)

points_growth = imodal.Utilities.fill_area_uniform_density(imodal.Utilities.area_shape, aabb_source, points_density, shape=1.5*shape_source)
points_growth_mask = imodal.Utilities.area_shape(points_growth, shape=1.*shape_source)
points_growth = points_growth[points_growth_mask]

# Initial normal frames for the growth module.
rot_growth = torch.stack([imodal.Utilities.rot2d(0.)]*points_growth.shape[0], axis=0)

Plot points of the growth module.

plt.plot(shape_source[:, 0].numpy(), shape_source[:, 1].numpy(), color='black')
plt.plot(points_growth[:, 0].numpy(), points_growth[:, 1].numpy(), 'o', color='blue')
plt.axis('equal')
plt.show()

Create the deformation model with a combination of 3 modules : implicit module of order 1 (growth model), implicit module of order 0 (small corrections) and a global translation.

Create and initialize the global translation module.

coeff_global = 15.

global_translation = imodal.DeformationModules.GlobalTranslation(2, coeff=coeff_global)

Create and initialize the growth module.

nu = 0.1
coeff_growth = 0.01
scale_growth = 100.

C = torch.empty(points_growth.shape[0], 2, 1)

growth = imodal.DeformationModules.ImplicitModule1(2, points_growth.shape[0], scale_growth, C, coeff=coeff_growth, nu=nu, gd=(points_growth, rot_growth))

Create and initialize local translations module.

coeff_small = 100.
nu = 0.1
scale_small = 10.
points_small = shape_source.clone()
small_scale_translations = imodal.DeformationModules.ImplicitModule0(2, points_small.shape[0], scale_small, coeff=coeff_small, nu=nu, gd=points_small)

Define our growth model tensor.

# The polynomial model for our growth model tensor.
def pol(pos, a, b, c, d):
    return a + b*pos[:, 1]  + c*pos[:, 1]**2 + d*pos[:, 1]**3


# Callback called when evaluating the model to compute the growth model tensor from parameters.
def callback_compute_c(init_manifold, modules, parameters, deformables):
    abcd = parameters['abcd']['params'][0]
    a = abcd[0].unsqueeze(1)
    b = abcd[1].unsqueeze(1)
    c = abcd[2].unsqueeze(1)
    d = abcd[3].unsqueeze(1)
    modules[3].C = pol(init_manifold[3].gd[0], a, b, c, d).transpose(0, 1).unsqueeze(2)


# Initial parameters of our growth model tensor.
abcd = torch.zeros(4, 2)
abcd[0] = 0.1 * torch.ones(2)
abcd.requires_grad_()

Out:

tensor([[0.1000, 0.1000],
        [0.0000, 0.0000],
        [0.0000, 0.0000],
        [0.0000, 0.0000]], requires_grad=True)

Define deformables used by the registration model.

deformable_shape_source = imodal.Models.DeformablePoints(shape_source)
deformable_shape_target = imodal.Models.DeformablePoints(shape_target)
deformable_dots_source = imodal.Models.DeformablePoints(dots_source)
deformable_dots_target = imodal.Models.DeformablePoints(dots_target)

Define the registration model.

model = imodal.Models.RegistrationModel(
    [deformable_shape_source, deformable_dots_source],
    [global_translation, growth, small_scale_translations],
    [imodal.Attachment.VarifoldAttachment(2, [50., 300.]),
     imodal.Attachment.EuclideanPointwiseDistanceAttachment(100.)],
    lam=200., other_parameters={'abcd': {'params': [abcd]}},
    precompute_callback=callback_compute_c)

Fitting using Torch LBFGS optimizer.

shoot_solver = 'euler'
shoot_it = 10

costs = {}
fitter = imodal.Models.Fitter(model, optimizer='torch_lbfgs')
fitter.fit([deformable_shape_target, deformable_dots_target], 200, costs=costs, options={'shoot_solver': shoot_solver, 'shoot_it': shoot_it, 'line_search_fn': 'strong_wolfe'})

Out:

Starting optimization with method torch LBFGS, using solver euler with 10 iterations.
Initial cost={'deformation': 0.0, 'attach': 389514912.0}
1e-10
Evaluated model with costs=389514912.0
Evaluated model with costs=389017248.0784309
Evaluated model with costs=386759138.37632895
Evaluated model with costs=360678089.23446655
Evaluated model with costs=262140285.0810547
Evaluated model with costs=121890011.91992188
Evaluated model with costs=93003823.87695312
Evaluated model with costs=92206315.44921875
Evaluated model with costs=90143611.37304688
Evaluated model with costs=81597550.7890625
Evaluated model with costs=72592011.3046875
Evaluated model with costs=68825231.87109375
Evaluated model with costs=67356862.42578125
Evaluated model with costs=66452007.34765625
Evaluated model with costs=66087941.10546875
Evaluated model with costs=63723252.7265625
Evaluated model with costs=58526261.38671875
Evaluated model with costs=53242483.0859375
Evaluated model with costs=42919517.75
Evaluated model with costs=39089258.54296875
Evaluated model with costs=36472452.53515625
Evaluated model with costs=35885266.71875
Evaluated model with costs=35845736.10546875
Evaluated model with costs=35809788.92578125
Evaluated model with costs=35618638.89453125
Evaluated model with costs=33920214.6015625
Evaluated model with costs=33916584.4765625
Evaluated model with costs=33873124.41015625
================================================================================
Time: 43.663380257003155
Iteration: 0
Costs
deformation=52000.41015625
attach=33821124.0
Total cost=33873124.41015625
1e-10
Evaluated model with costs=33873124.41015625
Evaluated model with costs=73905388.76953125
Evaluated model with costs=33842895.48828125
Evaluated model with costs=32958255.27734375
Evaluated model with costs=32759891.1796875
Evaluated model with costs=32742820.671875
Evaluated model with costs=32739666.7734375
Evaluated model with costs=32733289.67578125
Evaluated model with costs=32714473.75
Evaluated model with costs=32669626.9609375
Evaluated model with costs=32562321.8046875
Evaluated model with costs=32323923.3515625
Evaluated model with costs=31863071.5703125
Evaluated model with costs=31161121.2578125
Evaluated model with costs=29925769.70703125
Evaluated model with costs=28874945.9765625
Evaluated model with costs=28547800.046875
Evaluated model with costs=28444326.9921875
Evaluated model with costs=28420440.59375
Evaluated model with costs=28417674.5546875
Evaluated model with costs=28414425.12890625
Evaluated model with costs=28414247.765625
================================================================================
Time: 82.46674059799989
Iteration: 1
Costs
deformation=61635.765625
attach=28352612.0
Total cost=28414247.765625
1e-10
Evaluated model with costs=28414247.765625
Evaluated model with costs=28414309.77734375
Evaluated model with costs=28414281.92578125
Evaluated model with costs=28414311.796875
Evaluated model with costs=28414247.765625
================================================================================
Time: 91.62942164700144
Iteration: 2
Costs
deformation=61635.765625
attach=28352612.0
Total cost=28414247.765625
================================================================================
Optimisation process exited with message: Convergence achieved.
Final cost=28414247.765625
Model evaluation count=55
Time elapsed = 91.6297703900018

Compute optimized deformation trajectory.

intermediates = {}
with torch.autograd.no_grad():
    deformed = model.compute_deformed(shoot_solver, shoot_it, intermediates=intermediates)
    deformed_shape = deformed[0][0]
    deformed_dots = deformed[1][0]
deformed_growth = intermediates['states'][-1][3].gd[0]
deformed_growth_rot = intermediates['states'][-1][3].gd[1]
global_translation_controls = [control[2] for control in intermediates['controls']]
growth_controls = [control[3] for control in intermediates['controls']]

Plot results.

plt.subplot(1, 3, 1)
plt.title("Source")
plt.plot(shape_source[:, 0].numpy(), shape_source[:, 1].numpy(), '-', color='black')
plt.plot(dots_source[:, 0].numpy(), dots_source[:, 1].numpy(), '.', color='black')
plt.axis('equal')

plt.subplot(1, 3, 2)
plt.title("Deformed source")
plt.plot(deformed_shape[:, 0], deformed_shape[:, 1], '-', color='blue')
plt.plot(deformed_dots[:, 0], deformed_dots[:, 1], '.', color='blue')
plt.axis('equal')

plt.subplot(1, 3, 3)
plt.title("Deformed source and target")
plt.plot(shape_target[:, 0].numpy(), shape_target[:, 1].numpy(), '-', color='red')
plt.plot(dots_target[:, 0].numpy(), dots_target[:, 1].numpy(), '.', color='red')
plt.plot(deformed_shape[:, 0], deformed_shape[:, 1], '-', color='blue')
plt.plot(deformed_dots[:, 0], deformed_dots[:, 1], '.', color='blue')
# plt.plot(deformed_growth[:, 0], deformed_growth[:, 1], '.', color='blue')
plt.axis('equal')
plt.show()

Evaluate estimated growth model tensor.

learned_abcd = abcd.detach()
learned_C = pol(model.init_manifold[3].gd[0].detach(),
                learned_abcd[0].unsqueeze(1),
                learned_abcd[1].unsqueeze(1),
                learned_abcd[2].unsqueeze(1),
                learned_abcd[3].unsqueeze(1)).transpose(0, 1).unsqueeze(2).detach()
print("Learned growth model tensor parameters:\n {}".format(learned_abcd))

Out:

Learned growth model tensor parameters:
 tensor([[ 1.0000e-01,  1.0000e-01],
        [-1.6238e-07, -5.1295e-06],
        [-8.3878e-06,  2.0106e-04],
        [ 9.6801e-06,  8.4670e-06]])

Plot estimated growth factor.

ax = plt.subplot()
plt.plot(shape_source[:, 0].numpy(), shape_source[:, 1].numpy(), '-')
imodal.Utilities.plot_C_ellipses(ax, points_growth, learned_C, R=deformed_growth_rot, scale=1.)
plt.axis(aabb_source.squared().totuple())
plt.axis('equal')
plt.show()

Recompute the learned deformation trajectory this time with the grid deformation to visualize growth.

# We extract the modules of the models and fill the right manifolds.
modules = imodal.DeformationModules.CompoundModule(copy.copy(model.modules))
modules.manifold.fill(model.init_manifold.clone())
silent_shape = copy.copy(modules[0])
silent_dots = copy.copy(modules[1])
global_translation = copy.copy(modules[2])
growth = copy.copy(modules[3])

# Define the deformation grid.
square_size = 1.
growth_grid_resolution = [math.floor(aabb_source.width/square_size),
                          math.floor(aabb_source.height/square_size)]
deformation_grid = imodal.DeformationModules.DeformationGrid(aabb_source, growth_grid_resolution)

# We construct the controls we will give will shooting.
controls = [[torch.tensor([]), torch.tensor([]), torch.tensor([]), global_translation_control, growth_control] for growth_control, global_translation_control in zip(growth_controls, global_translation_controls)]

# Reshoot.
intermediates_growth = {}
with torch.autograd.no_grad():
    imodal.HamiltonianDynamic.shoot(imodal.HamiltonianDynamic.Hamiltonian([silent_shape, silent_dots, deformation_grid, global_translation, growth]), shoot_solver, shoot_it, controls=controls, intermediates=intermediates_growth)

# Store final deformation.
shoot_deformed_shape = silent_shape.manifold.gd.detach()
shoot_deformed_dots = silent_dots.manifold.gd.detach()
shoot_deformed_grid = deformation_grid.togrid()

Plot the deformation grid.

ax = plt.subplot()
plt.plot(shape_source[:, 0].numpy(), shape_source[:, 1].numpy(), '--', color='black')
plt.plot(dots_source[:, 0].numpy(), dots_source[:, 1].numpy(), '.', color='black')
plt.plot(shape_target[:, 0].numpy(), shape_target[:, 1].numpy(), '.-', color='red')
plt.plot(dots_target[:, 0].numpy(), dots_target[:, 1].numpy(), '.', color='black')
plt.plot(shoot_deformed_shape[:, 0].numpy(), shoot_deformed_shape[:, 1].numpy())
plt.plot(shoot_deformed_dots[:, 0].numpy(), shoot_deformed_dots[:, 1].numpy(), '.')
imodal.Utilities.plot_grid(ax, shoot_deformed_grid[0], shoot_deformed_grid[1], color='xkcd:light blue', lw=0.4)
plt.axis('equal')

plt.show()

Perform curve registration using the previously learned growth factor

Redefine deformation modules.

global_translation = imodal.DeformationModules.GlobalTranslation(2, coeff=coeff_global)

growth = imodal.DeformationModules.ImplicitModule1(2, points_growth.shape[0], scale_growth, learned_C, coeff=coeff_growth, nu=nu, gd=(points_growth, rot_growth))

small_scale_translation = imodal.DeformationModules.ImplicitModule0(2, shape_source.shape[0], scale_small, coeff=coeff_small, nu=nu, gd=shape_source)

Redefine deformables and registration model.

deformable_shape_source = imodal.Models.DeformablePoints(shape_source)
deformable_shape_target = imodal.Models.DeformablePoints(shape_target)

refit_model = imodal.Models.RegistrationModel([deformable_shape_source],
                [global_translation, growth, small_scale_translation],
                [imodal.Attachment.VarifoldAttachment(2, [50., 300.])],
                lam=100)

Fitting using Torch LBFGS optimizer.

shoot_solver = 'euler'
shoot_it = 10

costs = {}
fitter = imodal.Models.Fitter(refit_model, optimizer='torch_lbfgs')
fitter.fit([deformable_shape_target], 200, costs=costs, options={'shoot_solver': shoot_solver, 'shoot_it': shoot_it, 'line_search_fn': 'strong_wolfe'})

Out:

Starting optimization with method torch LBFGS, using solver euler with 10 iterations.
Initial cost={'deformation': 0.0, 'attach': 57479712.0}
1e-10
Evaluated model with costs=57479712.0
Evaluated model with costs=57384212.08517962
Evaluated model with costs=56520040.51559448
Evaluated model with costs=47695985.265625
Evaluated model with costs=26498847.1484375
Evaluated model with costs=9639712.80859375
Evaluated model with costs=5433212.6015625
Evaluated model with costs=3449238.1640625
Evaluated model with costs=2436832.2890625
Evaluated model with costs=1252332.609375
Evaluated model with costs=1118726.3203125
Evaluated model with costs=976751.8046875
Evaluated model with costs=816844.0859375
Evaluated model with costs=795104.671875
Evaluated model with costs=794498.375
Evaluated model with costs=792804.671875
Evaluated model with costs=790014.484375
Evaluated model with costs=781801.4296875
Evaluated model with costs=764769.828125
Evaluated model with costs=735108.1484375
Evaluated model with costs=690100.7421875
Evaluated model with costs=590372.2578125
Evaluated model with costs=544375.3515625
Evaluated model with costs=538265.0078125
================================================================================
Time: 41.00408576600239
Iteration: 0
Costs
deformation=70210.3203125
attach=468054.6875
Total cost=538265.0078125
1e-10
Evaluated model with costs=538265.0078125
Evaluated model with costs=531507.5234375
Evaluated model with costs=528703.5546875
Evaluated model with costs=526320.6328125
Evaluated model with costs=525346.484375
Evaluated model with costs=524065.2109375
Evaluated model with costs=521056.4296875
Evaluated model with costs=516694.28125
Evaluated model with costs=522931.8984375
Evaluated model with costs=513068.4765625
Evaluated model with costs=510881.171875
Evaluated model with costs=509233.03125
Evaluated model with costs=507304.7890625
Evaluated model with costs=504117.5703125
Evaluated model with costs=500183.25
Evaluated model with costs=489300.3203125
Evaluated model with costs=16499890.04296875
Evaluated model with costs=467400.6171875
Evaluated model with costs=1742351.0703125
Evaluated model with costs=476170.125
Evaluated model with costs=467066.46875
Evaluated model with costs=466209.53125
Evaluated model with costs=462545.4140625
Evaluated model with costs=461350.015625
Evaluated model with costs=458412.609375
================================================================================
Time: 81.97763606700028
Iteration: 1
Costs
deformation=96601.671875
attach=361810.9375
Total cost=458412.609375
1e-10
Evaluated model with costs=458412.609375
Evaluated model with costs=455251.7578125
Evaluated model with costs=445662.515625
Evaluated model with costs=433257.34375
Evaluated model with costs=422789.46875
Evaluated model with costs=430500.359375
Evaluated model with costs=421001.6015625
Evaluated model with costs=418981.421875
Evaluated model with costs=418639.65625
Evaluated model with costs=417944.4765625
Evaluated model with costs=416821.0
Evaluated model with costs=415191.15625
Evaluated model with costs=413834.1484375
Evaluated model with costs=412825.3046875
Evaluated model with costs=411898.640625
Evaluated model with costs=409821.9453125
Evaluated model with costs=406385.1796875
Evaluated model with costs=389744.671875
Evaluated model with costs=393897.859375
Evaluated model with costs=387350.2109375
Evaluated model with costs=385137.890625
Evaluated model with costs=384301.0078125
Evaluated model with costs=384227.71875
================================================================================
Time: 118.38591409199944
Iteration: 2
Costs
deformation=96871.46875
attach=287356.25
Total cost=384227.71875
1e-10
Evaluated model with costs=384227.71875
Evaluated model with costs=384101.3828125
Evaluated model with costs=383643.6015625
Evaluated model with costs=382828.890625
Evaluated model with costs=442592.2421875
Evaluated model with costs=382480.265625
Evaluated model with costs=381886.65625
Evaluated model with costs=381745.875
Evaluated model with costs=381380.15625
Evaluated model with costs=381688.2109375
Evaluated model with costs=380951.8515625
Evaluated model with costs=380696.9921875
Evaluated model with costs=380520.75
Evaluated model with costs=380442.0703125
Evaluated model with costs=378804.796875
Evaluated model with costs=371151.546875
Evaluated model with costs=386909.640625
Evaluated model with costs=369114.4609375
Evaluated model with costs=363798.71875
Evaluated model with costs=361630.109375
Evaluated model with costs=361396.6640625
Evaluated model with costs=361365.2734375
Evaluated model with costs=361215.8046875
Evaluated model with costs=360878.9609375
Evaluated model with costs=360127.8828125
================================================================================
Time: 157.24498883400156
Iteration: 3
Costs
deformation=101274.7578125
attach=258853.125
Total cost=360127.8828125
1e-10
Evaluated model with costs=360127.8828125
Evaluated model with costs=358982.2734375
Evaluated model with costs=357429.8046875
Evaluated model with costs=354586.671875
Evaluated model with costs=352394.5546875
Evaluated model with costs=351348.3828125
Evaluated model with costs=349321.25
Evaluated model with costs=347591.015625
Evaluated model with costs=345267.484375
Evaluated model with costs=345146.90625
Evaluated model with costs=344610.046875
Evaluated model with costs=343939.28125
Evaluated model with costs=343318.9765625
Evaluated model with costs=342172.0234375
Evaluated model with costs=336057.828125
Evaluated model with costs=330556.96875
Evaluated model with costs=329407.953125
Evaluated model with costs=328609.71875
Evaluated model with costs=328430.359375
Evaluated model with costs=328132.6875
Evaluated model with costs=327968.0703125
Evaluated model with costs=327889.9609375
================================================================================
Time: 191.51706939500218
Iteration: 4
Costs
deformation=92889.9609375
attach=235000.0
Total cost=327889.9609375
1e-10
Evaluated model with costs=327889.9609375
Evaluated model with costs=327868.875
Evaluated model with costs=327809.5
Evaluated model with costs=327817.8984375
Evaluated model with costs=327819.234375
Evaluated model with costs=327831.3671875
Evaluated model with costs=327823.4921875
Evaluated model with costs=327809.5
Evaluated model with costs=327809.5
Evaluated model with costs=327809.5
================================================================================
Time: 207.0719389320002
Iteration: 5
Costs
deformation=92962.625
attach=234846.875
Total cost=327809.5
1e-10
Evaluated model with costs=327809.5
Evaluated model with costs=327817.8984375
Evaluated model with costs=327819.234375
Evaluated model with costs=327831.3671875
Evaluated model with costs=327823.4921875
Evaluated model with costs=327809.5
Evaluated model with costs=327809.5
Evaluated model with costs=327809.5
================================================================================
Time: 219.47133496899914
Iteration: 6
Costs
deformation=92962.625
attach=234846.875
Total cost=327809.5
================================================================================
Optimisation process exited with message: Convergence achieved.
Final cost=327809.5
Model evaluation count=137
Time elapsed = 219.4716922459993

Compute optimized deformation trajectory.

intermediates = {}
with torch.autograd.no_grad():
    deformed = refit_model.compute_deformed(shoot_solver, shoot_it, intermediates=intermediates)
    deformed_shape = deformed[0][0]
deformed_growth = intermediates['states'][-1][2].gd[0]
deformed_growth_rot = intermediates['states'][-1][2].gd[1]
global_translation_controls = [control[1] for control in intermediates['controls']]
growth_controls = [control[2] for control in intermediates['controls']]
translation_controls = [control[3] for control in intermediates['controls']]

Plot results.

plt.subplot(1, 3, 1)
plt.title("Source")
plt.plot(shape_source[:, 0].numpy(), shape_source[:, 1].numpy(), '-', color='black')
plt.axis('equal')

plt.subplot(1, 3, 2)
plt.title("Deformed source")
plt.plot(deformed_shape[:, 0], deformed_shape[:, 1], '-', color='blue')
plt.axis('equal')

plt.subplot(1, 3, 3)
plt.title("Deformed source and target")
plt.plot(shape_target[:, 0].numpy(), shape_target[:, 1].numpy(), '-', color='red')
plt.plot(deformed_shape[:, 0], deformed_shape[:, 1], '-', color='blue')
plt.axis('equal')

plt.show()

Recompute the learned deformation trajectory this time with the grid deformation to visualize growth.

modules = imodal.DeformationModules.CompoundModule(copy.copy(refit_model.modules))
modules.manifold.fill(refit_model.init_manifold)

square_size = 1
grid_resolution = [math.floor(aabb_source.width/square_size),
                   math.floor(aabb_source.height/square_size)]
deformation_grid = imodal.DeformationModules.DeformationGrid(aabb_source, growth_grid_resolution)

controls = [control[1:] for control in intermediates['controls']]

deformable_shape = imodal.Models.DeformablePoints(shape_source)
deformable_shape.silent_module.manifold.cotan = refit_model.init_manifold[0].cotan
deformable_grid = imodal.Models.DeformableGrid(aabb_source, grid_resolution)

intermediates = {}
with torch.autograd.no_grad():
    imodal.Models.deformables_compute_deformed([deformable_shape, deformable_grid], modules[1:], shoot_solver, shoot_it, intermediates=intermediates, controls=controls)

Plot the growth trajectory.

indices = [0, 3, 7, 10]

fig = plt.figure(figsize=[5.*len(indices), 5.])
for i, index in enumerate(indices):
    state = intermediates['states'][index]

    ax = plt.subplot(1, len(indices), i + 1)
    deformable_grid.silent_module.manifold.fill_gd(state[1].gd)
    grid_x, grid_y = deformable_grid.silent_module.togrid()
    imodal.Utilities.plot_grid(ax, grid_x, grid_y, color='xkcd:light blue', lw=0.4)

    plt.plot(shape_source[:, 0].numpy(), shape_source[:, 1].numpy(), color='black')
    plt.plot(shape_target[:, 0].numpy(), shape_target[:, 1].numpy(), color='red')
    plt.plot(state[0].gd[:, 0].numpy(), state[0].gd[:, 1].numpy())

    plt.axis('equal')
    plt.axis('off')

fig.tight_layout()
plt.show()