Surface registration

Example of a diffeomorphic matching of surfaces using varifolds metrics: We perform an LDDMM matching of two meshes using the geodesic shooting algorithm.

Define our dataset

Standard imports

import os
import time

import imageio
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from mpl_toolkits.mplot3d import Axes3D
from torch.autograd import grad

from pykeops.torch import Kernel, kernel_product
from pykeops.torch.kernel_product.formula import *

# torch type and device
use_cuda = torch.cuda.is_available()
torchdeviceId = torch.device('cuda:0') if use_cuda else 'cpu'
torchdtype = torch.float32

# PyKeOps counterpart
KeOpsdeviceId = torchdeviceId.index  # id of Gpu device (in case Gpu is  used)
KeOpsdtype = torchdtype.__str__().split('.')[1]  # 'float32'

Import data file, one of :

  • “hippos.pt” : original data (6611 vertices),

  • “hippos_red.pt” : reduced size (1654 vertices),

  • “hippos_reduc.pt” : further reduced (662 vertices),

  • “hippos_reduc_reduc.pt” : further reduced (68 vertices)

if use_cuda:
    datafile = 'data/hippos.pt'
else:
    datafile = 'data/hippos_reduc_reduc.pt'

Define the kernels

Define Gaussian kernel \((K(x,y)b)_i = \sum_j \exp(-\gamma\|x_i-y_j\|^2)b_j\)

def GaussKernel(sigma):
    def K(x, y, b):
        params = {
            'id': Kernel('gaussian(x,y)'),
            'gamma': 1 / (sigma * sigma),
            'backend': 'auto'
        }
        return kernel_product(params, x, y, b)
    return K

Define “Gaussian-CauchyBinet” kernel \((K(x,y,u,v)b)_i = \sum_j \exp(-\gamma\|x_i-y_j\|^2) \langle u_i,v_j\rangle^2 b_j\)

def GaussLinKernel(sigma):
    def K(x, y, u, v, b):
        params = {
            'id': Kernel('gaussian(x,y) * linear(u,v)**2'),
            'gamma': (1 / (sigma * sigma), None),
            'backend': 'auto'
        }
        return kernel_product(params, (x, u), (y, v), b)
    return K

Custom ODE solver, for ODE systems which are defined on tuples

def RalstonIntegrator():
    def f(ODESystem, x0, nt, deltat=1.0):
        x = tuple(map(lambda x: x.clone(), x0))
        dt = deltat / nt
        l = [x]
        for i in range(nt):
            xdot = ODESystem(*x)
            xi = tuple(map(lambda x, xdot: x + (2 * dt / 3) * xdot, x, xdot))
            xdoti = ODESystem(*xi)
            x = tuple(map(lambda x, xdot, xdoti: x + (.25 * dt) * (xdot + 3 * xdoti), x, xdot, xdoti))
            l.append(x)
        return l

    return f

LDDMM implementation

Deformations: diffeomorphism

Hamiltonian system

def Hamiltonian(K):
    def H(p, q):
        return .5 * (p * K(q, q, p)).sum()
    return H


def HamiltonianSystem(K):
    H = Hamiltonian(K)
    def HS(p, q):
        Gp, Gq = grad(H(p, q), (p, q), create_graph=True)
        return -Gq, Gp
    return HS

Shooting approach

def Shooting(p0, q0, K, nt=10, Integrator=RalstonIntegrator()):
    return Integrator(HamiltonianSystem(K), (p0, q0), nt)


def Flow(x0, p0, q0, K, deltat=1.0, Integrator=RalstonIntegrator()):
    HS = HamiltonianSystem(K)
    def FlowEq(x, p, q):
        return (K(x, q, p),) + HS(p, q)
    return Integrator(FlowEq, (x0, p0, q0), deltat)[0]


def LDDMMloss(K, dataloss, gamma=0):
    def loss(p0, q0):
        p,q = Shooting(p0, q0, K)[-1]
        return gamma * Hamiltonian(K)(p0, q0) + dataloss(q)
    return loss

Data attachment term

Varifold data attachment loss for surfaces

# VT: vertices coordinates of target surface,
# FS,FT : Face connectivity of source and target surfaces
# K kernel
def lossVarifoldSurf(FS, VT, FT, K):
    def CompCLNn(F, V):
        V0, V1, V2 = V.index_select(0, F[:, 0]), V.index_select(0, F[:, 1]), V.index_select(0, F[:, 2])
        C, N = .5 * (V0 + V1 + V2), .5 * torch.cross(V1 - V0, V2 - V0)
        L = (N ** 2).sum(dim=1)[:, None].sqrt()
        return C, L, N / L

    CT, LT, NTn = CompCLNn(FT, VT)
    cst = (LT * K(CT, CT, NTn, NTn, LT)).sum()

    def loss(VS):
        CS, LS, NSn = CompCLNn(FS, VS)
        return cst + (LS * K(CS, CS, NSn, NSn, LS)).sum() - 2 * (LS * K(CS, CT, NSn, NTn, LT)).sum()

    return loss

Registration

Load the dataset and plot it

VS, FS, VT, FT = torch.load(datafile)
q0 = VS.clone().detach().to(dtype=torchdtype, device=torchdeviceId).requires_grad_(True)
VT = VT.clone().detach().to(dtype=torchdtype, device=torchdeviceId)
FS = FS.clone().detach().to(dtype=torch.long, device=torchdeviceId)
FT = FT.clone().detach().to(dtype=torch.long, device=torchdeviceId)
sigma = torch.tensor([20], dtype=torchdtype, device=torchdeviceId)

fig = plt.figure()
ax = Axes3D(fig)
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.plot_trisurf(q0.detach().cpu().numpy()[:, 0],
                q0.detach().cpu().numpy()[:, 1],
                q0.detach().cpu().numpy()[:, 2],
                triangles=FS.detach().cpu().numpy(),
                color=(0, 0, 0, 0),  edgecolor=(1, 0, 0, .08), linewidth=1)
ax.plot_trisurf(VT.detach().cpu().numpy()[:, 0],
                VT.detach().cpu().numpy()[:, 1],
                VT.detach().cpu().numpy()[:, 2],
                triangles=FT.detach().cpu().numpy(),
                color=(0, 0, 0, 0),  edgecolor=(0, 0, 1, .3),  linewidth=1)
blue_proxy = plt.Rectangle((0, 0), 1, 1, fc="b")
red_proxy = plt.Rectangle((0, 0), 1, 1, fc=(1, 0, 0, .5))
ax.legend([red_proxy,  blue_proxy], ['source', 'target'])
ax.set_title('Data')
plt.show()
../../_images/sphx_glr_plot_LDDMM_Surface_001.png

Out:

/home/bcharlier/tmp/libkeops/pykeops/tutorials/surface_registration/plot_LDDMM_Surface.py:208: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

Define data attachment and LDDMM functional

dataloss = lossVarifoldSurf(FS, VT, FT, GaussLinKernel(sigma=sigma))
Kv = GaussKernel(sigma=sigma)
loss = LDDMMloss(Kv, dataloss)

Out:

Compiling libKeOpstorche9c3dca64d in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorche9c3dca64d:
       formula: Sum_Reduction(((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * Pow((X_1|Y_1),2)) * B_0),0)
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); X_1 = Vi(2,3); Y_0 = Vj(3,3); Y_1 = Vj(4,3); B_0 = Vj(5,1);
       dtype  : float32
... Done.

Perform optimization

# initialize momentum vectors
p0 = torch.zeros(q0.shape, dtype=torchdtype, device=torchdeviceId, requires_grad=True)

optimizer = torch.optim.LBFGS([p0], max_eval=10, max_iter=10)
print('performing optimization...')
start = time.time()

def closure():
    optimizer.zero_grad()
    L = loss(p0, q0)
    print('loss', L.detach().cpu().numpy())
    L.backward()
    return L

for i in range(10):
    print('it ', i, ': ', end='')
    optimizer.step(closure)

print('Optimization (L-BFGS) time: ', round(time.time() - start, 2), ' seconds')

Out:

performing optimization...
it  0 : Compiling libKeOpstorchf44721a1c0 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorchf44721a1c0:
       formula: Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0)
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3);
       dtype  : float32
... Done.
Compiling libKeOpstorch288995f887 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch288995f887:
       formula: Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(1,3,0), Var(4,3,0), Var(5,3,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0);
       dtype  : float32
... Done.
Compiling libKeOpstorch3975a99230 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch3975a99230:
       formula: Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(2,3,1), Var(4,3,0), Var(5,3,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0);
       dtype  : float32
... Done.
Compiling libKeOpstorchece7967efa in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorchece7967efa:
       formula: Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(3,3,1), Var(4,3,0), Var(5,3,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0);
       dtype  : float32
... Done.
loss 91181.375
Compiling libKeOpstorch6f2b8d15fe in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch6f2b8d15fe:
       formula: Grad_WithSavedForward(Sum_Reduction(((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * Pow((X_1|Y_1),2)) * B_0),0), Var(1,3,0), Var(6,1,0), Var(7,1,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); X_1 = Vi(2,3); Y_0 = Vj(3,3); Y_1 = Vj(4,3); B_0 = Vj(5,1); Var(6,1,0); Var(7,1,0);
       dtype  : float32
... Done.
Compiling libKeOpstorchfd77c86ace in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorchfd77c86ace:
       formula: Grad_WithSavedForward(Sum_Reduction(((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * Pow((X_1|Y_1),2)) * B_0),0), Var(2,3,0), Var(6,1,0), Var(7,1,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); X_1 = Vi(2,3); Y_0 = Vj(3,3); Y_1 = Vj(4,3); B_0 = Vj(5,1); Var(6,1,0); Var(7,1,0);
       dtype  : float32
... Done.
Compiling libKeOpstorch889eb4a97f in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch889eb4a97f:
       formula: Grad_WithSavedForward(Sum_Reduction(((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * Pow((X_1|Y_1),2)) * B_0),0), Var(3,3,1), Var(6,1,0), Var(7,1,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); X_1 = Vi(2,3); Y_0 = Vj(3,3); Y_1 = Vj(4,3); B_0 = Vj(5,1); Var(6,1,0); Var(7,1,0);
       dtype  : float32
... Done.
Compiling libKeOpstorch0bb3b9486e in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch0bb3b9486e:
       formula: Grad_WithSavedForward(Sum_Reduction(((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * Pow((X_1|Y_1),2)) * B_0),0), Var(4,3,1), Var(6,1,0), Var(7,1,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); X_1 = Vi(2,3); Y_0 = Vj(3,3); Y_1 = Vj(4,3); B_0 = Vj(5,1); Var(6,1,0); Var(7,1,0);
       dtype  : float32
... Done.
Compiling libKeOpstorchf76f81171d in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorchf76f81171d:
       formula: Grad_WithSavedForward(Sum_Reduction(((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * Pow((X_1|Y_1),2)) * B_0),0), Var(5,1,1), Var(6,1,0), Var(7,1,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); X_1 = Vi(2,3); Y_0 = Vj(3,3); Y_1 = Vj(4,3); B_0 = Vj(5,1); Var(6,1,0); Var(7,1,0);
       dtype  : float32
... Done.
Compiling libKeOpstorchf259d43dec in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorchf259d43dec:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(3,3,1), Var(4,3,0), Var(5,3,0)), Var(1,3,0), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorch6da4d7fb62 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch6da4d7fb62:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(3,3,1), Var(4,3,0), Var(5,3,0)), Var(2,3,1), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorch4fbee33586 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch4fbee33586:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(3,3,1), Var(4,3,0), Var(5,3,0)), Var(3,3,1), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorch6ea557227e in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch6ea557227e:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(3,3,1), Var(4,3,0), Var(5,3,0)), Var(4,3,0), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorch3bc77f9d21 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch3bc77f9d21:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(2,3,1), Var(4,3,0), Var(5,3,0)), Var(1,3,0), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorcha9c613a16c in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorcha9c613a16c:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(2,3,1), Var(4,3,0), Var(5,3,0)), Var(2,3,1), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorch707d8c6227 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch707d8c6227:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(2,3,1), Var(4,3,0), Var(5,3,0)), Var(3,3,1), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorch9b9d9f027b in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch9b9d9f027b:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(2,3,1), Var(4,3,0), Var(5,3,0)), Var(4,3,0), Var(6,3,1), Var(7,3,1))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,1); Var(7,3,1);
       dtype  : float32
... Done.
Compiling libKeOpstorch4abd9d9140 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch4abd9d9140:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(1,3,0), Var(4,3,0), Var(5,3,0)), Var(1,3,0), Var(6,3,0), Var(7,3,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,0); Var(7,3,0);
       dtype  : float32
... Done.
Compiling libKeOpstorche7f48d2fd8 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorche7f48d2fd8:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(1,3,0), Var(4,3,0), Var(5,3,0)), Var(2,3,1), Var(6,3,0), Var(7,3,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,0); Var(7,3,0);
       dtype  : float32
... Done.
Compiling libKeOpstorche2dec47aa3 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorche2dec47aa3:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(1,3,0), Var(4,3,0), Var(5,3,0)), Var(3,3,1), Var(6,3,0), Var(7,3,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,0); Var(7,3,0);
       dtype  : float32
... Done.
Compiling libKeOpstorch6b3079c7a1 in /home/bcharlier/tmp/libkeops/pykeops/common/../build//build-libKeOpstorch6b3079c7a1:
       formula: Grad_WithSavedForward(Grad_WithSavedForward(Sum_Reduction((Exp( -(WeightedSqDist(G_0,X_0,Y_0))) * B_0),0), Var(1,3,0), Var(4,3,0), Var(5,3,0)), Var(4,3,0), Var(6,3,0), Var(7,3,0))
       aliases: G_0 = Pm(0,1); X_0 = Vi(1,3); Y_0 = Vj(2,3); B_0 = Vj(3,3); Var(4,3,0); Var(5,3,0); Var(6,3,0); Var(7,3,0);
       dtype  : float32
... Done.
loss 84040.125
loss 26011.625
loss 19569.25
loss 16383.25
loss 14535.3125
loss 10051.25
loss 9112.625
loss 7602.5
loss 7533.25
it  1 : loss 7533.25
loss 7459.8125
loss 7357.5625
loss 7194.8125
loss 6773.75
loss 6232.0
loss 5672.375
loss 5336.1875
loss 5147.3125
loss 4827.3125
it  2 : loss 4827.3125
loss 4333.8125
loss 3663.6875
loss 3325.9375
loss 2899.9375
loss 2759.625
loss 2681.875
loss 2521.5625
loss 2276.75
loss 2153.4375
it  3 : loss 2153.4375
loss 2100.5
loss 2073.5
loss 2034.5
loss 1959.8125
loss 1834.25
loss 1787.25
loss 1693.125
loss 1652.5
loss 1607.125
it  4 : loss 1607.125
loss 1538.3125
loss 1466.125
loss 1377.9375
loss 1327.875
loss 1305.25
loss 1287.75
loss 1252.875
loss 1223.1875
loss 1202.75
it  5 : loss 1202.75
loss 1189.25
loss 1176.875
loss 1165.8125
loss 1149.4375
loss 1137.25
loss 1118.625
loss 1081.875
loss 1042.5625
loss 1015.0
it  6 : loss 1015.0
loss 1008.375
loss 1001.1875
loss 999.6875
loss 997.75
loss 994.375
loss 986.125
loss 973.75
loss 960.125
loss 953.1875
it  7 : loss 953.1875
loss 949.4375
loss 946.6875
loss 944.125
loss 941.6875
loss 934.8125
loss 923.625
loss 901.625
loss 875.4375
loss 854.875
it  8 : loss 854.875
loss 847.875
loss 841.125
loss 837.9375
loss 834.375
loss 831.4375
loss 829.0625
loss 824.5625
loss 814.5
loss 803.6875
it  9 : loss 803.6875
loss 787.6875
loss 778.0
loss 774.0625
loss 772.875
loss 771.875
loss 770.8125
loss 769.75
loss 768.25
loss 766.0
Optimization (L-BFGS) time:  591.11  seconds

Display output

The animated version of the deformation:

nt = 15
listpq = Shooting(p0, q0, Kv, nt=nt)

The code to generate the .gif:

VTnp, FTnp = VT.detach().cpu().numpy(), FT.detach().cpu().numpy()
q0np, FSnp = q0.detach().cpu().numpy(), FS.detach().cpu().numpy()


images = []
for t in range(nt):
    qnp = listpq[t][1].detach().cpu().numpy()

    # create Figure
    fig = Figure(figsize=(6, 5), dpi=100)
    # Link canvas to fig
    canvas = FigureCanvasAgg(fig)

    # make the plot
    ax = Axes3D(fig)
    ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    ax.plot_trisurf(qnp[:, 0],  qnp[:, 1],  qnp[:, 2],  triangles=FSnp, color=(1, 1, 0, .5), edgecolor=(1, 1, 1, .3),  linewidth=1)
    ax.plot_trisurf(VTnp[:, 0], VTnp[:, 1], VTnp[:, 2], triangles=FTnp, color=(0, 0, 0, 0),  edgecolor=(0, 0, 1, .3),  linewidth=1)

    yellow_proxy = plt.Rectangle((0, 0), 1, 1, fc="y")
    ax.legend([yellow_proxy, blue_proxy], ['deformed', 'target'])
    ax.set_title('LDDMM matching example, step ' + str(t))

    # draw it!
    canvas.draw()

    # save plot in a numpy array through buffer
    s, (width, height) = canvas.print_to_buffer()
    images.append(np.frombuffer(s, np.uint8).reshape((height, width, 4)))

save_folder = '../../../doc/_build/html/_images/'
os.makedirs(save_folder, exist_ok=True)
imageio.mimsave(save_folder + 'surface_matching.gif', images, duration=.5)

Total running time of the script: ( 10 minutes 43.289 seconds)

Gallery generated by Sphinx-Gallery