Gaussian processes for non-rigid id regis istration - Con - - PowerPoint PPT Presentation

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Gaussian processes for non-rigid id regis istration

• Con
• nnections to
• medical im

image an analysis is

Marcel Lüthi

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

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Aims of the talk

• Show how Analysis by Synthesis and

Gaussian processes lead to a family of methods for non-rigid registration

• Provide an understanding of many

common algorithms in terms of Gaussian processes

• Show how to derive new registration

approaches using GPMMs and MCMC

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Outline

• The registration problem
• Problem formulation
• Registration as analysis by synthesis

problem

• An algorithm using Gaussian process

priors

• Priors for registration
• Spline-based models, Radial basis

functions

• Multis-scale and Spatially-varying models
• Statistical deformation models
• Likelihood functions
• Landmark registration
• Image to image registration
• Surface to image registration
• Advancing registration

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Outline

• The registration problem
• Problem formulation
• Registration as analysis by synthesis

problem

• An algorithm using Gaussian process

priors

• Priors for registration
• Spline-based models
• Radial basis functions
• Statistical deformation models
• Likelihood functions
• Landmark registration
• Image to image registration
• Surface to image registration
• Some ideas where to go from here

𝑦

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Ω

The registration problem

5

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

𝑦

𝜒(𝑦) Ω

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

The registration problem

Variational formulation 𝜒∗ = arg min

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

Mapping 𝜒∗ is trade-off that

• makes the images look similar (for similarity measure 𝐸)
• matches the prior assumptions (encoded by regularizer R)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Ω

The registration problem

7

Ω

Variational formulation 𝜒∗ = arg min

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

𝜒

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Ω

The registration problem

8

Ω

Variational formulation 𝜒∗ = arg min

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

𝜒

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Ω

The registration problem

9

Ω

Variational formulation 𝜒∗ = arg min

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

𝜒

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Representation of the mapping 𝜒

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Representation of the mapping 𝜒

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | 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 2016 | BASEL

Representation of the mapping 𝜒

Assumption: Images are rigidly aligned

• Mapping can be represented as a

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

• 𝜒 is parametric :

𝜒 𝜄 𝑦 = 𝑦 + 𝑣[𝜄](𝑦)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Representation of the mapping 𝜒

Mapping: 𝜒 𝜄 𝑦 = 𝑦 + 𝑣[𝜄](𝑦) Observation:

• Knowledge of 𝜄 and 𝐽𝑆 allows us to

synthesize target image 𝐽𝑈

• (at least up to intensity differences)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Registration as analysis by synthesis

Parameters 𝜄

Comparison: 𝑞 𝐽𝑈 𝜄, 𝐽𝑆) Update using 𝑞(𝜄|𝐽𝑈, 𝐽𝑆) Synthesis 𝜔(𝜄)

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

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Probabilistic formulation of registration

Mapping 𝜄∗ is trade-off that defines a mapping 𝜒[𝜄∗] which

• explains the data well (likelihood function)
• matches the prior assumptions (prior distribution)

MAP solution 𝜄∗ = arg max

𝜄

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

𝜄

𝑞 𝜒 𝜄 𝑞(𝐽𝑈|𝐽𝑆 ∘ 𝜒[𝜄]) Using Bayes rule: 𝑄 𝜒[𝜄]|𝐽𝑈, 𝐽𝑆 =

𝑄 𝐽𝑈|𝜒[𝜄],𝐽𝑆 𝑄 𝜒[𝜄],𝐽𝑆 𝑄 𝐽𝑈

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Registration problem

𝜒∗ = arg max

𝜄

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

𝜄

ln 𝑞(𝜒 𝜄 ) + ln(𝑞 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄] ) = arg min

𝜄 − ln 𝑞 𝜒 𝜄

− ln(𝑞 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄] )

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

The registration problem

Take home message: Registration is model fitting!!! Probabilistic formulation 𝜄∗ = arg min

𝜄 − ln 𝑞 𝐽𝑈 𝐽𝑆 ∘ 𝜒 𝜄

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

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

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Gaussian processes

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

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | 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 GP using only the first 𝑠 components of its KL-Expansion 𝑣 = 𝜈 + ෍

𝑗=1 𝑠

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

Parametric representation of Gaussian process

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Registration problem

𝜒∗ = arg min

𝜄 − ln 𝑞 𝜒 𝜄

− ln 𝑞(𝐽𝑈|𝐽𝑆 ∘ 𝜒[𝜄]) = arg min

𝜄 −ln 1

Z exp(− 1 2 𝜄 2) − ln(𝑞 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄] = arg min

𝜄

−ln 1 𝑎 + 1 2 𝜄 2 − ln(𝑞 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄] = arg min

𝜄

1 2 𝜄 2− ln(𝑞 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄]

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Summary: registration problem

• Variational and probabilistic formulation are closely related
• Prior can be seen as regularizer
• Likelihood term is an image similarity
• For a low-rank Gaussian process prior, the problem becomes parametric since

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

𝑗=1 𝑠

𝜄𝑗 𝜇𝑗 𝜚𝑗(𝑦)

• Can be optimized using gradient-descent schemes.
• All the regularization assumptions are encoded in the eigenfunctions 𝜚𝑗

arg min

𝜄

1 2 𝜄 2+ ln(𝑞 𝐽𝑈 𝐽𝑆 ∘ 𝜒[𝜄] )

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Outline

• The registration problem
• Problem formulation
• Registration as analysis by synthesis

problem

• An algorithm using Gaussian process

priors

• Priors for registration
• Spline-based models, Radial basis

functions

• Multis-scale and Spatially-varying models
• Statistical deformation models
• Likelihood functions
• Landmark registration
• Image to image registration
• Surface to image registration
• Advancing registration

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Why are priors interesting?

𝜄

𝜄∗ = arg max

𝜄

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

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Why are priors interesting?

𝜄∗ = arg max

𝜄

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

𝜄

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

ex

Defining a Gaussian process

A Gaussian process 𝐻𝑄 𝜈, 𝑙 is completely specified by a mean function 𝜈 and covariance function (or kernel) 𝑙.

• 𝜈: Ω → ℝ𝑒 defines how the average deformation looks like
• 𝑙: Ω × Ω → ℝ𝑒×𝑒 defines how it can deviate from the mean
• Must be positive semi-definite

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

The mean function

• Usual assumption:

𝜈 𝑦 = 𝜈1(𝑦) ⋮ 𝜈𝑒(𝑦) = ⋮

• The reference shape is an

average shape.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

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

𝑡 = 1, 𝜏 = 3 𝑡 = 1, 𝜏 = 5 𝑡 = 2, 𝜏 = 3

Scalar-valued Gaussian kernel

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

𝑙 𝑦, 𝑦′ = 𝑙(1)(𝑦, 𝑦′) ⋯ ⋮ ⋱ ⋮ ⋯ 𝑙(𝑒)(𝑦, 𝑦′)

• 𝑙(1), … , 𝑙(𝑒): 𝒴 × 𝒴 → ℝ are scalar-valued kernels
• 𝑙 ∶ 𝒴 × 𝒴 → ℝ𝑒×𝑒 becomes a matrix valued kernel.

Assumption: Each dimension is modelled independently.

• the output-dimensions are uncorrelated.

Diagonal kernel

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

𝑙 𝑦, 𝑦′ = s1exp − 𝑦 − 𝑦′ 2 𝜏1

2

s2exp − 𝑦 − 𝑦′ 2 𝜏2

2

A model for smooth 2D deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

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

A model for smooth deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

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

A model for smooth deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

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

A model for smooth deformations

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Matern class of kernels

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

𝜉

𝐿𝑤( 2𝜉 𝑦 − 𝑦′ 𝜍 )

• Γ is the Γ function, 𝑙𝜉 the modified Bessel function and 𝜉, 𝜍 are parameters
• The derivatives are 𝜉 − 1 times differentiable
• Special cases:
• 𝜉 =

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

)

• 𝜉 =

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

) exp(−

3 𝑦−𝑦′ 𝜍

)

• 𝜉 → ∞ Gaussian kernel

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Thin-plate splines

• Minimize the bending energy of a metal sheet

𝑆 𝜒 = ෍

𝑙=1 𝑒

න

Ω

𝛼𝑈𝛼𝜒𝑙 𝑦

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.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Elastic body splines

• Mechanical model of an elastic body or material
• Solution to the following PDE

𝜈𝛼2𝑣 𝑦 + 𝜈 + 𝜇 𝛼 𝛼 ⋅ 𝑣 𝑦 = 𝑑|𝑦|

• Corresponding (matrix-valued) covariance function (may not be positive definite)

𝑙 𝑦, 𝑦′ = 12 1 − 𝜉 − 1 |𝑦|2𝐽 − 3𝑦𝑦𝑈 where 𝜉 = 𝜇 2(𝜇 + 𝜈)

Kohlrausch, Jan, Karl Rohr, and H. Siegfried Stiehl. "A new class of elastic body splines for nonrigid registration of medical images." Journal of Mathematical Imaging and Vision 23.3 (2005): 253-280.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

B-Splines

• Use B-Spline basis function

𝑙 𝑦, 𝑧 = ෍

𝑙∈ℤ𝑒

𝛾 𝑡𝑦 − 𝑙 𝛾 𝑡𝑧 − 𝑙 Where 𝑡 is a scaling constant

• 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 2016 | 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 2016 | BASEL

Summary: Priors

• Leads to formulation of many standard transformation models in terms of Gaussian

process

• Improves understanding of methods
• Let’s us switch between “priors”
• Purely conceptual formulation
• No algorithms
• Can sample and visualize deformations
• Invaluable to check assumptions

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Outline

• The registration problem
• Problem formulation
• Registration as analysis by synthesis

problem

• An algorithm using Gaussian process

priors

• Priors for registration
• Spline-based models, Radial basis

functions

• Multis-scale and Spatially-varying models
• Statistical deformation models
• Likelihood functions
• Landmark registration
• Image to image registration
• Surface to image registration
• Advancing registration

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Landmark likelihood

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

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

𝑞 𝑚1

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

= ෑ

𝑗

𝑂 𝜒 𝜄 𝑚𝑆 , 𝜏2

𝑚1

𝑆

𝑚1

𝑈

𝑚𝑗

𝑆

𝑚𝑗

𝑈

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Given:

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

𝑆, ෤

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

• Find posterior distribution

𝑣 | 𝑚1

𝑆, … , 𝑚𝑜 𝑆, ෤

𝑣1, … , ෤ 𝑣𝑛

෤ 𝑣2

Landmark likelihood and GP Regression

෤ 𝑣1 𝑚𝑆

1

𝑚𝑆

𝑜

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | 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 2016 | BASEL

Landmark registration

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Landmark likelihood: Some remarks

• Classical problems 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 2016 | BASEL

Image to image registration

• What is a good synthesis function?

Simple choice: Use the warped reference image! 𝐽𝑆 ∘ ℎ𝜄

−1

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Image likelihood (single point)

• Probabilistic model: 𝐽𝑈 ℎ𝜄 𝑦

= 𝐽𝑆 𝑦 + 𝜗, 𝜗 ∼ 𝑂 0, 𝜏2 , 𝑦 ∈ ΩR

• Likelihood for a single point 𝑦:

𝑞 𝐽𝑈 ℎ𝜄(𝑦) 𝜄, 𝐽𝑆, 𝑦 ∼ 𝑂 𝐽𝑆 𝑦 , 𝜏2

𝐽𝑆: Ω → ℝ 𝐽𝑈: Ω → ℝ 𝐽𝑆(x) ℎ𝜄 𝐽𝑈(ℎ𝜄(x))

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Image likelihood (full image)

• Assuming that noise is independent at each point:

𝑞 𝐽𝑈 ∘ ℎ𝜄 Ԧ 𝜄, 𝐽𝑆 ∼ ෑ

𝑦∈I𝑆

𝑂 𝐽𝑆 𝑦 , 𝜏2 𝑞 𝐽𝑈 ∘ ℎ𝜄 Ԧ 𝜄, 𝐽𝑆 = 1 𝑎 ෑ

𝑦∈𝐽𝑆

exp(− 𝐽𝑆 𝑦 − 𝐽𝑈 ℎ𝜄 𝑦

2

𝜏2 )

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

The sum of squared distance metric

ln 𝑞 𝐽𝑈 ∘ ℎ𝜄 Ԧ 𝜄, 𝐽𝑆 = ln 1 𝑎 ෑ

𝑦∈𝐽𝑆

exp − 𝐽𝑆 𝑦 − 𝐽𝑈 ℎ𝜄 𝑦

2

𝜏2 = −𝑎1 ෍

𝑦∈𝐽𝑆

𝐽𝑆 𝑦 − 𝐽𝑈 ℎ𝜄 𝑦

2

= 𝐸𝑇𝑇𝐸[𝐽𝑆, 𝐽𝑈, ℎ𝜄]

• The sum of squared differences implements an independence assumption
• The parameter 𝜏2 becomes a weighting constant.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Likelihood from other metrics

• We can use any standard image metric 𝐸 𝐽𝑆, 𝐽𝑈, ℎ𝜄 to define a likelihood function:

𝑞 𝐽𝑈 𝜄, 𝐽𝑆 = 1 𝑎 exp −(𝐸 𝐽𝑆, 𝐽𝑈, ℎ𝜄 )

• Examples:
• Normalized cross correlations
• Mutual information
• …
• Makes it possible to reformulate any standard registration problem into the analysis by

synthesis framework.

• Special case of collective likelihood

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

What about surface registration?

Reference (surface): Γ𝑆 Γ𝑆 Γ𝑈 ℎ𝜄 Target (surface): Γ𝑈

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | 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 2016 | BASEL

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

Surface registration as image 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 2016 | BASEL

𝜍1 𝑦1

Active shape models (surface to image registration)

• ASMs model each profile as a normal

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

• Single profile point 𝑦𝑗:

𝑞 I𝑈(ℎ𝜄(𝑦𝑗))|𝜄, 𝑦𝑗 = 𝑂(𝜈𝑗, Σ𝑗)

• Image likelihood:

𝑞 I𝑈(ℎ𝜄(𝑦))|𝜄, Γ

𝑆 = ෑ 𝑗

𝑂(𝜈𝑗, Σ𝑗)

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Summary: Likelihood functions

• Synthesis function is often just a warp of a reference
• Works well if modality and dimensionality is the same
• Leads to very simple systems
• We get probabilistic interpretations of some standard metrics
• Makes assumptions more clear
• If we do not want full interpretation any metric can be turned into a likelihood function

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Outline

• The registration problem
• Problem formulation
• Registration as analysis by synthesis

problem

• An algorithm using Gaussian process

priors

• Priors for registration
• Spline-based models
• Radial basis functions
• Statistical deformation models
• Likelihood functions
• Landmark registration
• Image to image registration
• Surface to image registration
• Advancing registration
• More expressive priors
• Hybrid registration
• ASM using the Metropolis Hastings

algorithm

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

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

(scaling)

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

(or relationship)

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

(and relationship)

• 6. 𝑙 𝑦, 𝑦′ = 𝑙(𝜚 𝑦 , 𝜚 𝑦′ )

(domain warp)

Constructing s.p.d. kernels

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Multi-scale kernels

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

𝑗=0 𝑜

෍

𝑙∈ℤ𝑒

𝛾 2−𝑗𝑦 − 𝑙 𝛾 2−𝑗𝑧 − 𝑙 Opfer, Roland. "Multiscale kernels." Advances in computational mathematics 25.4 (2006): 357-380.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Multi-scale kernel

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | 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 2016 | BASEL

Anisotropic priors

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | 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 2016 | BASEL

Spatially-varying priors

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

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

The posterior 𝑣 | ෤ 𝑦1 , … , ෤ 𝑦𝑛, ෤ 𝑣1, … , ෤ 𝑣𝑛 is a Gaussian process

𝐻𝑄 𝜈𝑞, 𝑙𝑞

Its parameters are known analytically.

Landmark registration using Gaussian processes

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Landmark registration

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Hybrid registration

• We can now combine landmark registration with intensity:

1. Compute a posterior model using landmarks 𝐻𝑄(𝜈, 𝑙) 2. Use 𝐻𝑄(𝜈𝑞, 𝑙𝑞) as prior for registration with any image likelihood you like

• Example of Bayesian inference:

𝑞 𝑣 → 𝑞 𝑣 𝑀𝑆, 𝑀𝑈 → 𝑞 𝑣 𝑀𝑆, 𝑀𝑈, 𝐽𝑆, 𝐽𝑈

• Elegant solution to hybrid registration

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 2016 | BASEL

• Frees us from “tyranny of differentiability”
• Easy to integrate contours
• Let’s us model effects such as outliers, artifacts, … in principled ways
• Makes it possible to integrate results of bottom up proposals (landmark detectors)
• Let’s us reason about uncertainty of a solution

Draw a sample 𝑦′ from 𝑅(𝑦′|𝑦) Propose With probability 𝛽 = min

𝑄 𝒚′ 𝑄 𝒚 , 1

accept 𝒚′ as new sample Verify

Use Metropolis-Hastings for registration

67

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Call to arms

• Our MCMC scheme was designed for

really difficult problems

• 3D => 2D
• Complex illumination
• No scale
• Uncontrolled environment
• Let’s start together to tackle the

complicated problems in medical image analysis.

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2016 | BASEL

Call to arms

• Our MCMC scheme was designed for

really difficult problems

• 3D => 2D
• Complex illumination
• No scale
• Uncontrolled environment
• Let’s start together to tackle the

complicated problems in medical image analysis