# Computational Anatomy and PDE

## Outline

1. Background

2. EPDiff equation

3. New gradient flow PDE

# Background

## Part 1

MR camera

CT scanner

Radiology or medical imaging is in the middle of a digital or computer revolution that is affecting all of science and its various fields, including medicine

– Nick Bryan, Radiological Society of North America

• Better understanding of disease
• Tool to facilitate diagnosis

What is distance between images?

## Riemannian geometry: distance in curved space

Exemple: great circles (2d)

Exemple: Einstein's space-time (4d)

Exemple: space of diffeomorphisms (∞-dim)

\mathrm{Diff}(\Omega) = \{ \varphi\in C^{\infty}(\Omega,\Omega) \mid \varphi \text{ has smooth inverse}\}
$\mathrm{Diff}(\Omega) = \{ \varphi\in C^{\infty}(\Omega,\Omega) \mid \varphi \text{ has smooth inverse}\}$

## What is computational anatomy?

• Use Riemannian geometry on inf-dim manifold of diffeomorphisms to solve registration problems

(a) Template

(b) Warp

(c) Transformed

(d) Target

# EPDiff equation

## Diffeomorphic image registration

• Image domain $$\Omega$$
• Template image $$I_0: \Omega \to \mathbb{R}_+$$
• Target image $$I_1:\Omega\to \mathbb{R}_+$$
\displaystyle\min_{\varphi\in\mathrm{Diff}(\Omega)}\lVert I_0\circ\varphi^{-1}-I_1 \rVert^2
$\displaystyle\min_{\varphi\in\mathrm{Diff}(\Omega)}\lVert I_0\circ\varphi^{-1}-I_1 \rVert^2$
\varphi
$\varphi$

Ill-posed with noise!

## Large deformation diffeomorphic metric matching (LDDMM)

\displaystyle\min_{\varphi\in\mathrm{Diff}(\Omega)}\lVert I_0\circ\varphi^{-1}-I_1 \rVert^2 + \frac{1}{\sigma}\mathrm{dist}(\mathrm{id},\varphi)^2
$\displaystyle\min_{\varphi\in\mathrm{Diff}(\Omega)}\lVert I_0\circ\varphi^{-1}-I_1 \rVert^2 + \frac{1}{\sigma}\mathrm{dist}(\mathrm{id},\varphi)^2$
\displaystyle\mathrm{dist}(\mathrm{id},\varphi)^2 = \int_0^1 \lVert (1-\alpha \Delta)^k v(t) \rVert^2_{L^2} dt
$\displaystyle\mathrm{dist}(\mathrm{id},\varphi)^2 = \int_0^1 \lVert (1-\alpha \Delta)^k v(t) \rVert^2_{L^2} dt$
v(t) := \dot\varphi(t)\circ\varphi(t)^{-1} \quad \varphi(0)=\mathrm{id}, \varphi(1)=\varphi
$v(t) := \dot\varphi(t)\circ\varphi(t)^{-1} \quad \varphi(0)=\mathrm{id}, \varphi(1)=\varphi$
\displaystyle\min_{v\in C^\infty([0,1]\times\Omega,\mathbb{R}^n)}\lVert I_0\circ\varphi(1)^{-1}-I_1 \rVert^2 + \frac{1}{\sigma}\int_0^1 \lVert v(t,\cdot) \rVert^2_\mathcal{A} dt
$\displaystyle\min_{v\in C^\infty([0,1]\times\Omega,\mathbb{R}^n)}\lVert I_0\circ\varphi(1)^{-1}-I_1 \rVert^2 + \frac{1}{\sigma}\int_0^1 \lVert v(t,\cdot) \rVert^2_\mathcal{A} dt$
\frac{\partial}{\partial t}\varphi(t,\cdot) = v(t,\varphi(t,\cdot)), \quad \varphi(0,x) = x
$\frac{\partial}{\partial t}\varphi(t,\cdot) = v(t,\varphi(t,\cdot)), \quad \varphi(0,x) = x$

$$\varphi(1)$$ related to $$v$$ by

\lVert u \rVert_\mathcal{A} = \lVert \mathcal{A} u \rVert_{L^2} = \lVert (1-\alpha \Delta)^k u \rVert_{L^2}
$\lVert u \rVert_\mathcal{A} = \lVert \mathcal{A} u \rVert_{L^2} = \lVert (1-\alpha \Delta)^k u \rVert_{L^2}$

Norm $$|| \cdot||_\mathcal{A}$$ given by

## $$\Longrightarrow$$ EPDiff equation

• Velocity $$v:[0,1]\times\Omega\to\mathbb{R}^n$$
• Momentum $$m:[0,1]\times\Omega\to\mathbb{R}^n$$
\displaystyle\frac{\partial m}{\partial t} + \nabla m\cdot v + (\nabla v)^\top \cdot m + m(\nabla\cdot v) = 0
$\displaystyle\frac{\partial m}{\partial t} + \nabla m\cdot v + (\nabla v)^\top \cdot m + m(\nabla\cdot v) = 0$
\displaystyle m = (1-\alpha\Delta)^k v
$\displaystyle m = (1-\alpha\Delta)^k v$

## Numerical approach for LDDMM

• Steepest descent over time-dependent vector fields

• EPDiff IVP solver using particle methods (solitons)

• Use shooting over initial velocity $$v_0$$ to minimize
\displaystyle\lVert I_0\circ \varphi(1)^{-1}-I_1\rVert^2 + \frac{1}{\sigma}\lVert v_0 \rVert^2_\mathcal{A}
$\displaystyle\lVert I_0\circ \varphi(1)^{-1}-I_1\rVert^2 + \frac{1}{\sigma}\lVert v_0 \rVert^2_\mathcal{A}$

Result: excellent! (but computationally expensive)

## Open problems

• Develop FE discretization for EPDiff IVP
• Develop FE discretization for LDDMM
• Adaptive FE
• Numerical convergence analysis

*images by Matteo Molteni

# New gradient flow PDE

## Deformation gradient flows

Right-invariant metric on $$\mathrm{Diff}(\Omega)$$ defined by $$\mathcal A$$

Minimization problem

\displaystyle\min_{\varphi\in \mathrm{Diff}(\Omega)} \lVert(I_0\circ\varphi^{-1}-I_1\rVert^2
$\displaystyle\min_{\varphi\in \mathrm{Diff}(\Omega)} \lVert(I_0\circ\varphi^{-1}-I_1\rVert^2$
\displaystyle + \frac{1}{\sigma} d_{\text{FR}}(dx,\varphi_*dx)^2
$\displaystyle + \frac{1}{\sigma} d_{\text{FR}}(dx,\varphi_*dx)^2$

Origin of expensiveness: no formula for $$d_\mathrm{Diff}$$

\displaystyle + \frac{1}{\sigma} d_{\mathrm{Diff}}(\mathrm{id},\varphi)^2
$\displaystyle + \frac{1}{\sigma} d_{\mathrm{Diff}}(\mathrm{id},\varphi)^2$

Regularization penalizing change of volume

$$\mathrm{Diff}(\Omega)$$ acts on $$Q = C^{\infty}(\Omega)\times \mathrm{Dens}(\Omega)$$

$$E(\varphi) = F_{q_1}(\varphi\cdot q_0)$$ where $$q_i = (I_i,dx)$$

\dot \varphi = - \nabla_\mathcal{A} E(\varphi)
$\dot \varphi = - \nabla_\mathcal{A} E(\varphi)$
\mathrm{Diff}(\Omega)
$\mathrm{Diff}(\Omega)$
Q
$Q$
\mathrm{id}
$\mathrm{id}$
q_0
$q_0$
\varphi(t)
$\varphi(t)$
q_1
$q_1$
q(t) = \varphi(t)\cdot q_0
$q(t) = \varphi(t)\cdot q_0$
Q = C^{\infty}(\Omega)\times \mathrm{Dens}(\Omega)
$Q = C^{\infty}(\Omega)\times \mathrm{Dens}(\Omega)$
q_i = (I_i,dx)
$q_i = (I_i,dx)$

Gradient flow is

\dot \varphi = - v(\varphi\cdot q_0)\circ \varphi
$\dot \varphi = - v(\varphi\cdot q_0)\circ \varphi$
(1-\alpha\Delta)^s v = \frac{1}{\sigma}\rho \nabla \rho + \nabla(I-I_\infty)
$(1-\alpha\Delta)^s v = \frac{1}{\sigma}\rho \nabla \rho + \nabla(I-I_\infty)$
\rho = \det(\varphi^{-1})
$\rho = \det(\varphi^{-1})$
I = I_0\circ\varphi^{-1}
$I = I_0\circ\varphi^{-1}$

Lie-Euler method

\rho_k = \det(\varphi_k^{-1}) \quad I_k = I_0\circ\varphi_k^{-1}
$\rho_k = \det(\varphi_k^{-1}) \quad I_k = I_0\circ\varphi_k^{-1}$
\varphi_{k+1} = \exp(-h v_k)\circ\varphi_k
$\varphi_{k+1} = \exp(-h v_k)\circ\varphi_k$
(1-\alpha\Delta)^s v_k = \frac{1}{\sigma}\rho_k \nabla \rho_k + \nabla(I_k-I_\infty)
$(1-\alpha\Delta)^s v_k = \frac{1}{\sigma}\rho_k \nabla \rho_k + \nabla(I_k-I_\infty)$

## Open problem

• Develop moving mesh FE discretization for gradient flow PDE

Course project

slides available at slides.com/kmodin/cmm-kristineberg

By Klas Modin

# Computational Anatomy and PDE

Guest seminar in the Nordic Graduate Course on Computational Mathematical Modeling at Kristineberg on the west-coast of Sweden.

• 1,577