Computational Anatomy and PDE

Klas Modin


  1. Background
  2. EPDiff equation
  3. New gradient flow PDE


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}\}
Diff(Ω)={φC(Ω,Ω)φ has smooth inverse}\mathrm{Diff}(\Omega) = \{ \varphi\in C^{\infty}(\Omega,\Omega) \mid \varphi \text{ has smooth inverse}\}

Image dist = Riemannian dist between \(\varphi\) and \(\mathrm{id}\)

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

History of computational anatomy

EPDiff equation

Part 2

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
minφDiff(Ω)I0φ1I12\displaystyle\min_{\varphi\in\mathrm{Diff}(\Omega)}\lVert I_0\circ\varphi^{-1}-I_1 \rVert^2

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
minφDiff(Ω)I0φ1I12+1σdist(id,φ)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
dist(id,φ)2=01(1αΔ)kv(t)L22dt\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):=φ˙(t)φ(t)1φ(0)=id,φ(1)=φ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
minvC([0,1]×Ω,Rn)I0φ(1)1I12+1σ01v(t,)A2dt\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
tφ(t,)=v(t,φ(t,)),φ(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}
uA=AuL2=(1αΔ)kuL2\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

Euler-Lagrange equations

\(\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
mt+mv+(v)m+m(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
m=(1αΔ)kv\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}
I0φ(1)1I12+1σv0A2\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

Part 3

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
minφDiff(Ω)(I0φ1I12\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
+1σdFR(dx,φ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
+1σdDiff(id,φ)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)
φ˙=AE(φ)\dot \varphi = - \nabla_\mathcal{A} E(\varphi)
q(t) = \varphi(t)\cdot q_0
q(t)=φ(t)q0q(t) = \varphi(t)\cdot q_0
Q = C^{\infty}(\Omega)\times \mathrm{Dens}(\Omega)
Q=C(Ω)×Dens(Ω)Q = C^{\infty}(\Omega)\times \mathrm{Dens}(\Omega)
q_i = (I_i,dx)
qi=(Ii,dx)q_i = (I_i,dx)

Gradient flow is







\dot \varphi = - v(\varphi\cdot q_0)\circ \varphi
φ˙=v(φq0)φ\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αΔ)sv=1σρρ+(II)(1-\alpha\Delta)^s v = \frac{1}{\sigma}\rho \nabla \rho + \nabla(I-I_\infty)
\rho = \det(\varphi^{-1})
ρ=det(φ1)\rho = \det(\varphi^{-1})
I = I_0\circ\varphi^{-1}
I=I0φ1I = I_0\circ\varphi^{-1}

Lie-Euler method






\rho_k = \det(\varphi_k^{-1}) \quad I_k = I_0\circ\varphi_k^{-1}
ρk=det(φk1)Ik=I0φk1\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
φk+1=exp(hvk)φ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αΔ)svk=1σρkρk+(IkI)(1-\alpha\Delta)^s v_k = \frac{1}{\sigma}\rho_k \nabla \rho_k + \nabla(I_k-I_\infty)

Toy example (2D)

Open problem

  • Develop moving mesh FE discretization for gradient flow PDE

Course project

slides available at

Computational Anatomy and PDE

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