Non-rigid Registration Marcel Lthi Graphics and Vision Research - - PowerPoint PPT Presentation

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Non-rigid Registration

Marcel Lüthi

Graphics and Vision Research Group Department of Mathematics and Computer Science University of Basel

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Outline

• Non-rigid registration: The basic formulation
• Exercise: Parametric registration in Scalismo
• Advanced Priors
• Likelihood functions
• Exercise: ASMs in Scalismo
• Optimization

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Ω

The registration problem

3

Reference: 𝐽𝑆: Ω → ℝ Target: 𝐽𝑈: Ω → ℝ 𝜒: Ω → Ω

𝑦

𝜒(𝑦) Ω

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Why is it important?

• Do automatic measurements
• Compare shapes
• Statistics
• Build statistical models
• Transfer labels and annotations
• Atlas based segmentation

Ω 𝜒: Ω → Ω Ω

Maybe the most important problem in computer vision and medical image analysis

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Registration as analysis by synthesis

Parameters 𝜄

Comparison: 𝑞 𝐽𝑈 𝜄, 𝐽𝑆) Update using 𝑞(𝜄|𝐽𝑈, 𝐽𝑆) Synthesis 𝜒[𝜄]

Prior 𝜒[𝜄] ∼ 𝑞(𝜄) 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄]

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

The registration problem

Mapping 𝜒[𝜄∗] is trade-off that

• how well does the mapping explain the target image

(likelihood function)

• matches the prior assumptions (prior distribution)

Probabilistic formulation 𝜄∗ = arg max

𝜄

𝑞 𝜄 𝐽𝑈, 𝐽𝑆 = arg max

𝜄

𝑞 𝜄 𝑞(𝐽𝑈|𝜄, 𝐽𝑆)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Ω

The registration problem

Ω

𝜒[𝜄]

𝜄∗ = arg max

𝜄

𝑞 𝜄 𝐽𝑈, 𝐽𝑆 = arg max

𝜄

𝑞 𝜄 𝑞(𝐽𝑈|𝜄, 𝐽𝑆)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Ω

The registration problem

Ω 𝜄∗ = arg max

𝜄

𝑞 𝜄 𝐽𝑈, 𝐽𝑆 = arg max

𝜄

𝑞 𝜄 𝑞(𝐽𝑈|𝜄, 𝐽𝑆)

𝜒[𝜄]

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Ω

The registration problem

Ω

𝜒[𝜄]

𝜄∗ = arg max

𝜄

𝑞 𝜄 𝐽𝑈, 𝐽𝑆 = arg max

𝜄

𝑞 𝜄 𝑞(𝐽𝑈|𝜄, 𝐽𝑆)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

The registration problem

Main questions:

• How do we represent the mapping?
• How do we define the prior?
• What is the likelihood function?
• How can we solve the optimization problem?

Probabilistic formulation 𝜒∗ = arg max

𝜒

𝑞 𝜒 𝐽𝑈, 𝐽𝑆 = arg max

𝜒

𝑞 𝜒 𝑞(𝐽𝑈|𝜒, 𝐽𝑆)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Representation of the mapping 𝜒

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Representation of the mapping 𝜒

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Representation of the mapping 𝜒

Assumption: Images are rigidly aligned

• Mapping can be represented as a

displacement vector field: 𝜒 𝑦 = 𝑦 + 𝑣 𝑦 𝑣 ∶ Ω → ℝ𝑒

𝑦 u(𝑦)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Representation of the mapping 𝜒

Assumption: Images are rigidly aligned

• Mapping can be represented as a

displacement vector field: 𝜒 𝑦 = 𝑦 + 𝑣 𝑦 𝑣 ∶ Ω → ℝ𝑒 Observation: Knowledge of 𝑣 and 𝐽𝑆 allows us to synthesize target image 𝐽𝑈

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Registration as analysis by synthesis

Parameters 𝜄

Comparison: 𝑞 𝐽𝑈 𝜄, 𝐽𝑆) Update using 𝑞(𝜄|𝐽𝑈, 𝐽𝑆) Synthesis 𝜒[𝜄]

Prior 𝜒[𝜄] ∼ 𝑞(𝜄) 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄]

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Priors

Define the Gaussian process 𝑣 ∼ 𝐻𝑄 𝜈, 𝑙 with mean function 𝜈: Ω → ℝ2 and covariance function 𝑙: Ω × Ω → ℝ2×2 .

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Zero mean: 𝜈 𝑦 = 0 Squared exponential covariance function (Gaussian kernel) 𝑙 𝑦, 𝑦′ = s1exp − 𝑦 − 𝑦′ 2 𝜏1

2

s2exp − 𝑦 − 𝑦′ 2 𝜏2

2

Example prior: Smooth 2D deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

𝑡1 = 𝑡2 small, 𝜏1 = 𝜏2 large

Example prior: Smooth 2D deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

𝑡1 = 𝑡2 small, 𝜏1 = 𝜏2 small

Example prior: Smooth 2D deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

𝑡1 = 𝑡2 large, 𝜏1 = 𝜏2 large

Example prior: Smooth 2D deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

• We have a finite, parametric representation of the process.
• We know the pdf for a deformation 𝑣

𝑞 𝑣[𝛽] = 𝑞 𝛽 = ෑ

𝑗=1 𝑠

1 2𝜌 exp(−𝛽𝑗

2/2) = 1

𝑎 exp(− 1 2 𝛽 2) Represent 𝐻𝑄(𝜈, 𝑙) using only the first 𝑠 components of its KL-Expansion 𝑣 = 𝜈 + ෍

𝑗=1 𝑠

𝛽𝑗 𝜇𝑗 𝜚𝑗, 𝛽𝑗 ∼ 𝑂(0, 1)

Parametric representation of Gaussian process

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Registration as analysis by synthesis

Parameters 𝜄

Comparison: 𝑞 𝐽𝑈 𝜄, 𝐽𝑆) Update using 𝑞(𝜄|𝐽𝑈, 𝐽𝑆) Synthesis 𝜒[𝜄]

Prior 𝜒[𝜄] ∼ 𝑞(𝜄) 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄]

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Likelihood function: Image registration

Assumptions:

• Corresponding points have the same image intensity (up to i.i.d. noise)

𝑦 𝜒[𝜄](𝑦) 𝐽𝑆 𝐽𝑈 𝑞 𝐽𝑈(𝜒[𝜄](𝑦)) 𝐽𝑆, 𝜄, 𝑦 ∼ 𝑂 𝐽𝑆 𝑦 , 𝜏2 Images are similar when the intensities match

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Likelihood function: Image registration

𝑦 𝜒[𝜄](𝑦) 𝐽𝑆 𝐽𝑈 𝑞 𝐽𝑈 𝐽𝑆, 𝜄 = ෑ

𝑦∈Ω

𝑞 𝐽𝑈(𝜒[𝜄](𝑦)) 𝐽𝑆, 𝜄, 𝑦 = ෑ

𝑦∈Ω

1 𝑎 exp − (𝐽𝑈 𝜒 𝑦 − 𝐽𝑆 𝑦 ) 2 𝜏2 Images are similar when the intensities match Assumptions:

• Corresponding points have the same image intensity (up to i.i.d. noise)

Image term outside mapping function. Makes problem really difficult

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Registration as analysis by synthesis

Parameters 𝜄

Comparison: 𝑞 𝐽𝑈 𝜄, 𝐽𝑆) Update using 𝑞(𝜄|𝐽𝑈, 𝐽𝑆) Synthesis 𝜒[𝜄]

Prior 𝜒[𝜄] ∼ 𝑞(𝜄) 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄]

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Registration problem

𝛽∗ = 𝜄∗ = arg max

𝜄

𝑞 𝜒 𝜄 𝑞(𝐽𝑈|𝜒[𝜄], 𝐽𝑆) = arg max

𝜄

1 𝑎1 exp − 1 2 𝜄 2 1 𝑎2 ෑ

𝑦

exp − 𝐽𝑈 𝜒 𝜄 𝑦 − 𝐽𝑆(𝑦))

2

𝜏2

• Parametric problem, since:

𝜒[𝜄](𝑦) = 𝑦 + 𝜈(𝑦) + ෍

𝑗=1 𝑠

𝜄𝑗 𝜇𝑗 𝜚𝑗(𝑦)

• Can be optimized using gradient descent

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Variational formulation

𝛽∗ = arg max

𝜄

1 𝑎1 exp − 1 2 𝜄 2 1 𝑎2 ෑ

𝑦

exp − 𝐽𝑈 𝜒 𝜄 𝑦 − 𝐽𝑆(𝑦))

2

𝜏2 arg max

𝜄

ln 1 𝑎1 exp − 1 2 𝜄 2 + ln 1 𝑎2 ෑ

𝑦

exp − 𝐽𝑈 𝜒 𝜄 𝑦 − 𝐽𝑆(𝑦))

2

𝜏2 = arg max

𝜄

ln 1 𝑎1 − 1 2 𝜄 2 + ln 1 𝑎2 − ෍

𝑦∈Ω

𝐽𝑈 𝜒 𝜄 𝑦 − 𝐽𝑆(𝑦))

2

𝜏2 = arg min

𝜄

෍

𝑦∈Ω

𝐽𝑈 𝜒 𝜄 𝑦 − 𝐽𝑆(𝑦))

2

𝜏2 + 𝜇 2 𝜄 2

Image metric Regularizer

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

The registration problem

Probabilistic formulation 𝜄∗ = arg min

𝜄 − ln 𝑞 𝐽𝑈 𝐽𝑆, 𝜒 𝜄

− ln 𝑞 𝜒 𝜄 Variational formulation 𝜄∗ = arg min

𝜄 𝐸 𝐽𝑈, 𝐽𝑆, 𝜒[𝜄] + 𝜇𝑆[𝜒 𝜄 ]

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Exercise: Registration in Scalismo

Type into the codepane: goto(“http://shapemodelling.cs.unibas.ch/exercises/Exercise14.html”) Scalismo 0.16: Check examples in https://github.com/unibas-gravis/pmm2018

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

A selection of Gaussian process priors

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Why are priors interesting?

𝜄

𝜄∗ = arg max

𝜄

𝑞 𝜒 𝜄 𝑞(𝐽𝑈|𝐽𝑆, 𝜒[𝜄])

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Why are priors interesting?

𝜄∗ = arg max

𝜄

𝑞 𝜒 𝜄 𝑞(𝐽𝑈|𝐽𝑆, 𝜒[𝜄])

𝜄

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Models of smooth deformations

• Typical assumption:
• Deformation field is smooth
• GP approach
• Choose smooth kernel functions

𝑙 𝑦, 𝑦′ = 𝑡 exp(− 𝑦 − 𝑦′ 2 𝜏2 )

• Regularization operators
• Penalize large derivatives

ℛ 𝑣 = 𝑆𝑣 2 = ෍

𝑗=0 𝑜

𝛽𝑗 𝐸𝑗𝑣 2

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Connection between regularizer and kernel

Discrete setting: Finite difference operators Regularizer: ℛ ො 𝑣 = ෡ 𝐸ො 𝑣

2

෡ 𝐸ො 𝑣 = −2 1 1 −2 1 1 −2 1 1 −2 1 1 −2 1 ⋱ 𝑣(𝑦1) 𝑣(𝑦2) 𝑣(𝑦3) 𝑣(𝑦4) 𝑣(𝑦5) ⋮

Steinke, Florian, and Bernhard Schölkopf. "Kernels, regularization and differential equations." Pattern Recognition 41.11 (2008): 3271-3286. 𝑣(𝑦1) 𝑣(𝑦𝑗) 𝑦1 𝑦𝑗

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Connection between regularizer and kernel

p ො u = exp(− 1 2 ℛ ො 𝑣 ) = exp(− 1 2 ෡ 𝐸ො 𝑣

2)

= exp(− 1 2 ෡ 𝐸ො 𝑣

𝑈(෡

𝐸ො 𝑣) = exp(− 1 2 ො 𝑣𝑈 ෡ 𝐸𝑈 ෡ 𝐸ො 𝑣 )

෡ 𝐸𝑈 ෡ 𝐸 = 𝐿−1

• D specifies the “inverse” covariance matrix
• Can compute K from:

෡ 𝐸𝑈 ෡ 𝐸𝐿 = 𝐽

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Green’s functions and covariance functions

ℛ 𝑣 = 𝑆𝑣 2 = ෍

𝑗=0 𝑜

𝛽𝑗 𝐸𝑗𝑣 2 Corresponding covariance function for GP is the Greens function G: 𝑆∗𝑆𝐻 𝑦, 𝑧 = 𝜀(𝑦 − 𝑧)

• We can define Gaussian processes, which mimic typical regularization operators.
• T. Poggio and F. Girosi; Networks for Approximation and Learning, Proceedings of the IEEE, 1990

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Example: Gaussian kernel

𝑙 𝑦, 𝑦′ = exp(− 𝑦 − 𝑦′ 2 𝜏2 )

ℛ 𝑣 =

𝑆𝑣 2 = ෍

𝑗=0 ∞ 𝜏2𝑗

𝑗! 2𝑗 𝐸𝑗𝑣 2

• Non-zero functions are penalized
• pushes functions to zero away from data

Yuille, A. and Grzywacz M. A mathematical analysis of the motion coherence theory. International Journal of Computer vision

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Example: Exponential kernel (1D case)

𝑙 𝑦, 𝑦′ = 1 2𝛽 exp(−𝛽 𝑦 − 𝑦′ ) ℛ 𝑣 = 𝑆𝑣 2 = 𝛽2𝑣 + 𝐸1𝑣 2

Rasmussen, Carl Edward, and Christopher KI Williams. Gaussian processes for machine learning. Vol. 1. Cambridge: MIT press, 2006.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Matérn class of kernels

𝑙 𝑦, 𝑦′ = 𝑡 21−𝜉 Γ 𝜉 2 2𝜉 𝑦 − 𝑦′ 𝜍

𝜉

𝐿𝑤( 2𝜉 𝑦 − 𝑦′ 𝜍 )

• Γ is the Γ function, 𝑙𝜉 the modified Bessel function and 𝜉, 𝜍 are parameters
• Process 𝑣~𝐻𝑄 0, 𝑙 is 𝜉 − 1 times m.s. differentiable
• Special cases:
• 𝜉 =

1 2 : 𝑙 𝑦, 𝑦′ = 𝑡 exp(− 𝑦−𝑦′ 𝜍

)

• 𝜉 =

3 2 : k x, x′ = 𝑡(1 + 3 𝑦−𝑦′ 𝜍

) exp(−

3 𝑦−𝑦′ 𝜍

)

• 𝜉 → ∞ Gaussian kernel

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Thin-plate splines

• Minimizes the bending energy of a metal sheet

𝑆 𝑣 = 𝛼𝑈𝛼𝑣

2

• Corresponding covariance function

𝑙 𝑦, 𝑦′ = 1 12 2 𝑦 − 𝑦′ 3 − 3𝑆( 𝑦 − 𝑦′ 2 + 𝑆3 where 𝑆 = max

𝑦,𝑦′∈Ω ‖𝑦 − 𝑦′‖ Rohr, Karl, et al. "Landmark-based elastic registration using approximating thin-plate splines." IEEE Transactions on medical imaging 20.6 (2001): 526-534.

Williams, Oliver and Fitzgibbon Andrew, “Gaussian process implicit surfaces”

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

B-Splines

• We can build a covariance function from B-Spline basis functions 𝛾

(𝑡 is a scaling constant) 𝑙 𝑦, 𝑧 = ෍

𝑙∈ℤ𝑒

𝛾 𝑡𝑦 − 𝑙 𝛾 𝑡𝑧 − 𝑙

• Corresponding deformation model often called “free form deformations”
• Rueckert, Daniel, et al. "Nonrigid registration using free-form deformations: application to breast MR images."

IEEE transactions on medical imaging 18.8 (1999): 712-721.

• Klein, Stefan, et al. "Elastix: a toolbox for intensity-based medical image registration." IEEE transactions on

medical imaging 29.1 (2010): 196-205.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Many standard models for registration can be formulated using Gaussian processes

• Yields probabilistic interpretation
• We can sample and visualize deformation fields
• Can use them as building blocks for more complicated priors

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

• 1. 𝑙 𝑦, 𝑦′ = 𝑔 𝑦 𝑔 𝑦’ 𝑈, 𝑔: 𝑌 → ℝ𝑒
• 2. 𝑙 𝑦, 𝑦′ = 𝛽𝑙1 𝑦, 𝑦′ , 𝛽 ∈ ℝ+

(scaling)

• 3. k 𝑦, 𝑦′ = 𝐶𝑈𝑙1 𝑦, 𝑦′ 𝐶, B ∈ ℝ𝑠×𝑒 (lifting)
• 4. 𝑙 𝑦, 𝑦′ = 𝑙1 𝑦, 𝑦′ + 𝑙2 𝑦, 𝑦′

(or relationship)

• 5. 𝑙 𝑦, 𝑦′ = 𝑙1 𝑦, 𝑦′ ⋅ 𝑙2(𝑦, 𝑦′)

(and relationship)

Constructing s.p.d. kernels

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Multi-scale kernels

Add kernels that act on different scales: 𝑙 𝑦, 𝑦′ = ෍

𝑗=0 𝑜

෍

𝑙∈ℤ𝑒

𝛾 2−𝑗𝑦 − 𝑙 𝛾 2−𝑗𝑦′ − 𝑙

• Wavelet like multiscale representation

Opfer, Roland. "Multiscale kernels." Advances in computational mathematics 25.4 (2006): 357-380.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Multi-scale kernel

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Anisotropic priors

Scale deformations differently in each direction k 𝑦, 𝑦′ = 𝑆𝑈 𝑡1 𝑡2 𝑙 𝑦, 𝑦′ 𝑡1 𝑡2 𝑆

• R is a rotation matrix
• 𝑙 is scalar valued
• 𝑡1, s2 scaling factors

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Anisotropic priors

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Spatially-varying priors

Use different models for different regions 𝑙 𝑦, 𝑦′ = 𝜓 𝑦 𝜓 𝑦′ 𝑙1 𝑦, 𝑦′ + 1 − 𝜓 𝑦 (1 − 𝜓 𝑦′ ) 𝑙2(𝑦, 𝑦′) χ 𝑦 = ቊ1 if 𝑦 ∈ thumb region

• therwise

𝜓 𝑦 = 1 𝜓 𝑦 = 0 Freiman, Moti, Stephan D. Voss, and Simon K. Warfield. "Demons registration with local affine adaptive regularization: application to registration of abdominal structures." Biomedical Imaging: From Nano to Macro, 2011 IEEE International Symposium on. IEEE, 2011.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Spatially-varying priors

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Statistical deformation models

Estimate mean and covariance function from data:

𝜈 𝑦 = 𝑣 𝑦 = 1 𝑜 ෍

𝑗−1 𝑜

𝑣𝑗 (𝑦) 𝑙𝑇𝑁 𝑦, 𝑦′ = 1 𝑜 − 1 ෍

𝑗 𝑜

(𝑣𝑗 𝑦 − 𝑣(𝑦)) 𝑣𝑗 𝑦′ − 𝑣(𝑦′)

𝑈

𝑣1 ∶ Ω → ℝ2 𝑣2 ∶ Ω → ℝ2 … 𝑣𝑜 ∶ Ω → ℝ2

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Example 5: Statistical deformation models

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

A selection of likelihood functions

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Landmark likelihood

For one landmark pair (𝑚𝑆, 𝑚𝑈): 𝑞 𝑚𝑈 𝜄, 𝑚𝑆 = 𝑂 𝜒 𝜄 𝑚𝑆 , 𝐽2𝑦2𝜏2 For many landmarks 𝑀 = ((𝑚𝑆

1, 𝑚𝑈 1), … , (𝑚𝑆 𝑜, 𝑚𝑈 𝑜))

𝑞 𝑚1

𝑈, … , 𝑚𝑜 𝑈 𝜄, 𝑚𝑆 1, … , 𝑚𝑆 𝑜

= ෑ

𝑗

𝑂 𝜒 𝜄 𝑚𝑆 , 𝐽2𝑦2𝜏2

𝑚1

𝑆

𝑚1

𝑈

𝑚𝑗

𝑆

𝑚𝑗

𝑈

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Landmark likelihood: Some remarks

• Classical problem in registration
• Needs either many landmark points or good structure of prior to achieve good results

Rohr, Karl, et al. "Landmark-based elastic registration using approximating thin-plate splines." IEEE Transactions on medical imaging 20.6 (2001): 526-534.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Given:

• Gaussian process: 𝑣 ∼ 𝐻𝑄(𝜈, 𝑙)
• Observations: {(𝑚𝑗

𝑆, ෤

𝑣𝑗), 𝑗 = 1 , … , 𝑜} Assume: ෤ 𝑣𝑗 = 𝑣 𝑚𝑗 + 𝜗 with 𝜗 ∼ 𝑂(0, 𝜏2𝐽2×2). Goal:

• Find posterior distribution

𝑣 | 𝑚1

𝑆, … , 𝑚𝑜 𝑆, ෤

𝑣1, … , ෤ 𝑣𝑜

𝑣𝑜

Landmark registration using GP Regression

𝑣1 𝑚𝑆

1

𝑚𝑆

𝑜

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

𝜈𝑞(𝑦) = 𝜈 𝑦 + 𝐿 𝑦, 𝑍 (𝐿 𝑍, 𝑍 + 𝜏2𝐽2𝑜×2𝑜 )−1 ෥ 𝒗 − 𝜈(𝑍) 𝑙𝑞 𝑦, 𝑦′ = 𝑙 𝑦, 𝑦′ − 𝐿 𝑦, 𝑍 (𝐿 𝑍, 𝑍 + 𝜏2𝐽2𝑜×2𝑜 )−1𝐿(𝑍, 𝑦′)

The posterior 𝑣 |𝑚1

𝑆, … , 𝑚𝑜 𝑆, ෤

𝑣1, … , 𝑜 is a Gaussian process

𝐻𝑄 𝜈𝑞, 𝑙𝑞

Its parameters are known analytically.

Gaussian process regression

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Landmark registration using GP Regression

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Hybrid registration

• Combine landmark registration with intensity:
• 1. Use Gaussian process regression to obtain posterior from 𝐻𝑄(𝜈, 𝑙) from landmarks
• 2. Use 𝐻𝑄(𝜈𝑞, 𝑙𝑞) as new prior model for registration
• Simple solution to otherwise difficult problem

Wörz, Stefan, and Karl Rohr. "Hybrid spline-based elastic image registration using analytic solutions of the navier equation." Bildverarbeitung für die Medizin 2007. Springer Berlin Heidelberg, 2007. 151-155. Lu, Huanxiang, Philippe C. Cattin, and Mauricio Reyes. "A hybrid multimodal non-rigid registration of MR images based on diffeomorphic demons." Engineering in Medicine and Biology Society (EMBC), 2010 Annual International Conference of the IEEE. IEEE, 2010.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Likelihood function: Image registration

Assumptions:

• Corresponding points have the same image intensity (up to i.i.d. noise)

𝑦 𝜒[𝜄](𝑦) 𝐽𝑆 𝐽𝑈 𝑞 𝐽𝑈(𝜒[𝜄](𝑦)) 𝐽𝑆, 𝜄, 𝑦 ∼ 𝑂 𝐽𝑆 𝑦 , 𝜏2 Images are similar when the intensities match

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Likelihood function: Image registration

𝑦 𝜒[𝜄](𝑦) 𝐽𝑆 𝐽𝑈 𝑞 𝐽𝑈 𝐽𝑆, 𝜄 = ෑ

𝑦∈Ω

𝑞 𝐽𝑈(𝜒[𝜄](𝑦)) 𝐽𝑆, 𝜄, 𝑦 = ෑ

𝑦∈Ω

1 𝑎 exp − (𝐽𝑈 𝜒 𝑦 − 𝐽𝑆 𝑦 ) 2 𝜏2 Images are similar when the intensities match Assumptions:

• Corresponding points have the same image intensity (up to i.i.d. noise)

Image term outside mapping function. Makes problem really difficult

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Image vs. Landmark registration

• Landmark registration is easy
• All components are Gaussian
• Closed form solution using Gaussian process regression
• Image registration is hard
• Image destroys Gaussian assumption
• Likelihood function is not Gaussian
• Problem with many local minima

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

What about surface registration?

Reference (surface): Γ𝑆 Γ𝑆 Γ𝑈 𝜒 Target (surface): Γ𝑈

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

A trick: Implicit definition of a surface

• Any surface Γ can be represented as

the zero level set of a level set function Φ Γ = 𝑦 Φ 𝑦 = 0}

• Popular choice is the signed distance

function defined as 𝐸Γ 𝑦 = CP

Γ(𝑦) − 𝑦

with CP

Γ(𝑦) = arg min 𝑦′∈Γ ‖𝑦 − 𝑦′‖ 𝐸Γ 𝑦 = −30 𝐸Γ 𝑦 = −15 𝐸Γ 𝑦 = 0

𝐸Γ 𝑦 = 15 𝐸Γ 𝑦 = 30 𝐸Γ 𝑦 = 45

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Reference 𝐸𝑆: Ω𝑆 → ℝ Target 𝐸𝑈 ∶ Ω𝑈 → ℝ

Likelihood function: Surface registration

• We define the distance functions and use image to image likelihoods

𝑞 𝐸𝑈 𝜒[𝜄](𝑦) 𝜄, 𝐸𝑆, 𝑦 ∼ 𝑂 𝐸𝑆𝜒[𝜄](𝑦), 𝜏2

• Most likely solution will map points on the zero level-sets to each other
• Noise parameter 𝜏2 has geometric interpretation (variance of distance between the mapped

points)

𝑦

𝜒[𝜄](𝑦)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

𝜍1 𝑦1

Likelihood function: Active shape models

• ASMs model each profile 𝜍(𝑦𝑗) as a normal

distribution 𝑞 𝜍(𝑦𝑗) = 𝑂(𝜈𝑗, Σ𝑗)

• Single profile point 𝑦𝑗:

𝑞 𝜍(𝜒[𝜄](𝑦𝑗))|𝜄, 𝑦𝑗 = 𝑂(𝜈𝑗, Σ𝑗)

• Image likelihood:

𝑞 𝜍(𝜒[𝜄](𝑦))|𝜄, Γ

𝑆 = ෑ 𝑗

𝑂(𝜈𝑗, Σ𝑗)

Extracts profile (feature) from image

Shape is well matched if environment around profile points is likeli under trained model.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Exercise: Active Shape Models in Scalismo

Type into the codepane: goto(“http://shapemodelling.cs.unibas.ch/exercises/Exercise16.html”) Scalismo 0.16: Check examples in https://github.com/unibas-gravis/pmm2018

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Summary: Modelling

• GPs provide probabilistic interpretation to classic registration models
• Can visualize assumptions
• New ways to combine priors to individual applications.
• Modelling of prior is separate from likelihood
• Any GP model can be combined with any likelihood function
• Flexible framework to tailor model and algorithm to

needs of applications

• No increase in complexity
• Modelling and model fitting are separated
• Gradient based methods are important, but not the only method

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Optimization

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

The optimization problem

• The final problem is a difficult
• ptimization problem
• Possibly many local minima
• Non-linearity due to image term
• Not possible to avoid it
• Flexible models makes things worse

𝜄∗ = arg max

𝜄

𝑞 𝜒 𝜄 𝑞(𝐽𝑈|𝐽𝑆 ∘ 𝜒[𝜄])

𝜄

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Local minima

• Rigid Transformation
• Minima due to structure of object

Possible approach: Multi-resolution

• Non-rigid Transformation
• Minima appear/dissappear when shape changes

Possible approach: Multi-scale models, regularization

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Multi resolution

Idea: Solve optimzation problem for a sequence of smoothed out objects.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Implementation

• Smooth the input shapes
• For images, achieved by Gaussian blurring

72 Almost no local minima No-details Many local minima All-details Initial registration Final registration

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Multi-scale / Regularization

Idea: Solve optimzation problem for a sequence of increasingly detailed deformations Only large, smooth deformations Large regularization value Allow detailed deformations Almost no regularization Initial registration Final registration

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Doing the registration

Strategies:

• Gradient-based registration
• Compute gradient and use local optimization methods
• Quasi-Newton schemes , SGD, …
• Gradient free registration
• Use global optimization method directly on cost function
• Examples: Simulated annealing, Particle Swarm, …
• ICP-based methods
• Assume correspondence and solve in each iteration analytic problem
• Examples: Non-rigid ICP, Active Shape models, CPD

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Model-fitting using Markov Chain Monte Carlo

• Can obtain full posterior distribution

using the Metropolis Hastings algorithm

• Needs only point-wise evaluation of

unnormlized posterior

• Leads to principled way to integrate

unreliable bottom up methods

• Automatically detected landmarks

MAP Solution 𝜒∗ = arg max

𝜒

𝑞 𝜒 𝐽𝑈, 𝐽𝑆

𝑞 𝜒 𝐽𝑈, 𝐽𝑆) = 𝑞 𝜒[𝛽] 𝑞 𝐽𝑈 𝐽𝑆, 𝜒[𝛽] 𝑂(𝐽𝑈)

More on this tomorrow!