Image Registration

* * * * *

Carlton CHU

Smooth Realign Normalise Segment DARTEL With slides by John Ashburner, Chloe Hutton and Jesper Andersson

Overview of SPM Analysis

fMRI time-series Design matrix Statistical Parametric Map Motion Correction Smoothing Spatial Normalisation General Linear Model Parameter Estimates Anatomical Reference

* * *

Contents

Preliminaries *

Smooth

* * * Rigid-Body and Affine Transformations Optimisation and Objective Functions Transformations and Interpolation Intra-Subject Registration Inter-Subject Registration

Smooth

Smoothing is done by convolution .

Each voxel after smoothing effectively becomes the result of applying a weighted region of interest (ROI).

Before convolution Convolved with a circle Convolved with a Gaussian

Image Registration

•

Registration

- i.e. Optimise the parameters that describe a spatial transformation between the source and reference (template) images •

Transformation

to the determined transformation parameters - i.e. Re-sample according

2D Affine Transforms

* Translations by t x and t y * * * x 1 y 1 = x 0 = y 0 + t x + t y Rotation around the origin by  * * * * * x 1 y 1 x 1 y 1 = cos(  ) x 0 = -sin(  ) x 0 Zooms by s x = s x x 0 = s y y 0 + sin(  ) y 0 + cos(  ) y 0 and s y radians * Shear * x 1 = x 0 * y 1 = y 0 + h y 0

2D Affine Transforms

* Translations by t x and t y * * * x 1 y 1 = 1 x 0 = 0 x 0 + 0 y 0 + 1 y 0 + t x + t y Rotation around the origin by  * * * * * x 1 y 1 x 1 y 1 = cos(  ) x 0 = -sin(  ) x 0 Zooms by s x + sin(  ) y 0 + cos(  ) y 0 and s y : = s x x 0 + 0 y 0 + 0 = 0 x 0 + s y y 0 + 0 + 0 + 0 radians * Shear * x 1 = 1 x 0 * y 1 = 0 x 0 + h y 0 + 0 + 1 y 0 + 0

3D Rigid-body Transformations

* * A 3D rigid body transform is defined by: * 3 translations - in X, Y & Z directions * 3 rotations - about X, Y & Z axes The order of the operations matters

1

   

0 0 0 0 1 0 0 0 0 1 0 X trans Y trans Zt rans

   

1 1

    

0 0 0

Translations

0 cos

Φ 

sin

Φ

0 0 sin

Φ

cos

Φ

0

Pitch about x axis

0 0 0

   

1



cos

Θ     

0 sin

Θ

0 0 1 0 0 sin

Θ

0 cos

Θ

0

Roll about y axis

0 0 0

    

1

    

cos

Ω

sin 0

Ω

0 sin

Ω

cos

Ω

0 0 0 0 1 0

Yaw about z axis

0 0 0

   

1

Voxel-to-world Transforms

* * Affine transform associated with each image * Maps from voxels (x=1..n

x , y=1..n

y , z=1..n

z ) to some world co-ordinate system. e.g., * Scanner co-ordinates - images from DICOM toolbox * T&T/MNI coordinates - spatially normalised Registering image B (source) to image A (target) will update B’s voxel-to-world mapping * Mapping from voxels in A to voxels in B is by * A-to-world using M A , then world-to-B using M B -1 * M B -1 M A

Optimisation

* * Optimisation involves finding some “best” parameters according to an “objective function”, which is either minimised or maximised The “objective function” is often related to a probability based on some model Objective function Most probable solution (global optimum) Local optimum Value of parameter Local optimum

Objective Functions

* * Intra-modal * Mean squared difference (minimise) * * Normalised cross correlation (maximise) Entropy of difference (minimise) Inter-modal (or intra-modal) * Mutual information (maximise) * * * Normalised mutual information (maximise) Entropy correlation coefficient (maximise) AIR cost function (minimise)

Transformation

* * Images are re-sampled. An example in 2D: for y 0 =1..n

y0 % loop over rows for x 0 =1..n

x0 x 1 y 1 = t x (x 0 ,y 0 ,q) % transform according to q = t y (x 0 ,y 0 ,q) if 1  x 1  % loop over pixels in row n x1 & 1  y 1  n y1 then % voxel in range f 1 (x 0 ,y 0 ) = f 0 (x 1 ,y 1 ) % assign re-sampled value end % voxel in range end % loop over pixels in row end % loop over rows What happens if x 1 and y 1 are not integers?

Simple Interpolation

* * Nearest neighbour * Take the value of the closest voxel Tri-linear * Just a weighted average of the neighbouring voxels * * * f 5 f 6 f 7 = f 1 = f 3 = f 5 x 2 x 2 y 2 + f 2 + f 4 + f 6 x 1 x 1 y 1

B-spline Interpolation

A continuous function is represented by a linear combination of basis functions 2D B-spline basis functions of degrees 0, 1, 2 and 3 B-splines are piecewise polynomials Nearest neighbour and trilinear interpolation are the same as B-spline interpolation with degrees 0 and 1.

* * *

Contents

Preliminaries Intra-Subject Registration *

Realign

* Mean-squared difference objective function * * Residual artifacts and distortion correction

Coregister

Inter-Subject Registration

Mean-squared Difference

* * Minimising mean-squared difference works for intra-modal registration (realignment) Simple relationship between intensities in one image, versus those in the other * Assumes normally distributed differences

* * * *

Residual Errors from aligned fMRI

Re-sampling can introduce interpolation errors * especially tri-linear interpolation Gaps between slices can cause aliasing artefacts Slices are not acquired simultaneously * rapid movements not accounted for by rigid body model Image artefacts may not move according to a rigid body model * image distortion * * image dropout Nyquist ghost * Functions of the estimated motion parameters can be modelled as confounds in subsequent analyses

Movement by Distortion Interaction of fMRI

* Subject disrupts B encode direction 0 field, rendering it inhomogeneous => distortions in phase * * Subject moves during EPI time series Distortions vary with subject orientation => shape varies

Movement by distortion interaction

Correcting for distortion changes using Unwarp

Estimate reference from mean of all scans.

Estimate movement parameters.

Estimate new distortion fields for each image: • estimate rate of change of field with respect to the current estimate of movement parameters in pitch and roll.

Unwarp time series.

 

B

0   +  

B

0  

Andersson et al, 2001

* * *

Contents

Preliminaries Intra-Subject Registration *

Realign

*

Coregister

* Mutual Information objective function Inter-Subject Registration

Inter-modal registration

• Match images from same subject but different modalities: – anatomical localisation of single subject activations – achieve more precise spatial normalisation of functional image using anatomical image.

Mutual Information

* * * Used for between-modality registration Derived from joint histograms MI= *  ab P(a,b) log 2 [P(a,b)/( P(a) P(b) )] Related to entropy: MI = -H(a,b) + H(a) + H(b) * Where H(a) =  a P(a) log 2 P(a) and H(a,b) =  a P(a,b) log 2 P(a,b)

* * *

Contents

Preliminaries Intra-Subject Registration Inter-Subject Registration *

Normalise

* Affine Registration * * Nonlinear Registration Regularisation *

Segment

Spatial Normalisation - Reasons

* * Inter-subject averaging * Increase sensitivity with more subjects * Fixed-effects analysis * Extrapolate findings to the population as a whole * Mixed-effects analysis Standard coordinate system * e.g., Talairach & Tournoux space

Spatial Normalisation - Procedure

* Minimise mean squared difference from template image(s) Affine registration Non-linear registration

T2 T1 Transm T1 305 T1 EPI PD Template Images PET PD T2 “Canonical” images SS PET A wider range of contrasts can be registered to a linear combination of template images.

PD Spatial normalisation can be weighted so that non brain voxels do not influence the result.

Similar weighting masks can be used for normalising lesioned brains.

Spatial Normalisation - Templates

Spatial Normalisation - Affine

* * The first part is a 12 parameter affine transform * 3 translations * * * 3 rotations 3 zooms 3 shears Fits overall shape and size * Algorithm simultaneously minimises * Mean-squared difference between template and source image * Squared distance between parameters and their expected values (regularisation)

Spatial Normalisation - Non-linear

Deformations consist of a linear combination of smooth basis functions These are the lowest frequencies of a 3D discrete cosine transform (DCT) Algorithm simultaneously minimises * Mean squared difference between template and source image * Squared distance between parameters and their known expectation

Spatial Normalisation - Overfitting

Without regularisation, the non-linear spatial normalisation can introduce unnecessary warps.

Template image Non-linear registration (  2 using regularisation.

= 302.7) (  2 Affine registration.

= 472.1) Non-linear registration (  without regularisation.

2 = 287.3)

* * *

Contents

Preliminaries Intra-Subject Registration Inter-Subject Registration *

Normalise

*

Segment

* Gaussian mixture model * * Intensity non-uniformity correction Deformed tissue probability maps

Segmentation

* Segmentation in SPM5 also estimates a spatial transformation that can be used for spatially normalising images.

* It uses a generative model , which involves: * * * Mixture of Gaussians (MOG) Bias Correction Component Warping (Non-linear Registration) Component

*

Gaussian Probability Density

If intensities are assumed to be Gaussian of mean m k and variance s 2 k , then the probability of a value y i is:

*

Non-Gaussian Probability Distribution

A non-Gaussian probability density function can be modelled by a Mixture of Gaussians (MOG): Mixing proportion - positive and sums to one

Belonging Probabilities

Belonging probabilities are assigned by normalising to one.

*

Mixing Proportions

The mixing proportion g k represents the prior probability of a voxel being drawn from class k irrespective of its intensity.

* So:

*

Non-Gaussian Intensity Distributions

Multiple Gaussians per tissue class allow non-Gaussian intensity distributions to be modelled.

* E.g. accounting for partial volume effects

*

Probability of Whole Dataset

If the voxels are assumed to be independent, then the probability of the whole image is the product of the probabilities of each voxel: * A maximum-likelihood solution can be found by minimising the negative log-probability:

*

Modelling a Bias Field

A bias field is included, such that the required scaling at voxel i , parameterised by b , is r i ( b ) .

* * Replace the means by m k / r i ( b ) Replace the variances by ( s k / r i ( b )) 2

Modelling a Bias Field

* After rearranging...

y

r ( b )

y

r ( b )

Tissue Probability Maps

* Tissue probability maps (TPMs) are used instead of the proportion of voxels in each Gaussian as the prior.

ICBM Tissue Probabilistic Atlases

Mazziotta and Arthur W. Toga.

. These tissue probability maps are kindly provided by the International Consortium for Brain Mapping, John C.

“Mixing Proportions”

* * Tissue probability maps for each class are included.

The probability of obtaining class k at voxel i , given weights g is then:

Deforming the Tissue Probability Maps

* * Tissue probability images are deformed according to parameters a .

The probability of obtaining class k at voxel i , given weights g and parameters a is then:

The Extended Model

* By combining the modified P(c i =k|  ) and P(y i |c i =k,  ) , the overall objective function ( E ) becomes:

The Objective Function

* * *

Optimisation

The “best” parameters are those that minimise this objective function.

Optimisation involves finding them.

Begin with starting estimates, and repeatedly change them so that the objective function decreases each time.

Steepest Descent

Start Optimum Alternate between optimising different groups of parameters

Schematic of optimisation

Repeat until convergence… Hold g , m , s 2 and a constant, and minimise E - Levenberg-Marquardt strategy, using dE/d b w.r.t. b and d 2 E/d b 2 Hold g , m , s 2 and b constant, and minimise E - Levenberg-Marquardt strategy, using dE/d a w.r.t. a and d 2 E/d a 2 end Hold a and b constant, and minimise E w.r.t. g , m and s 2 -Use an Expectation Maximisation (EM) strategy.

* * *

Levenberg-Marquardt Optimisation

LM optimisation is used for nonlinear registration ( a ) and bias correction ( b ).

Requires first and second derivatives of the objective function ( E ).

Parameters a and b are updated by * Increase l to improve stability (at expense of decreasing speed of convergence).

*

Expectation Maximisation is used to update

m

,

s 2

and

g For iteration (n) , alternate between: * E-step : Estimate belonging probabilities by: * M-step : Set  (n+1) to values that reduce:

* *

Regularisation

Some bias fields and warps are more probable (a priori) than others.

Encoded using Bayes rule (for a maximum a posteriori solution): * Prior probability distributions modelled by a multivariate normal distribution.

* Mean vector m a and m b * * * Covariance matrix S a and S b -log[P( a )] = ( a m a ) T S a -1 ( am a ) -log[P( b )] = ( b m b ) T S b -1 ( bm b ) + const + const

Initial Affine Registration

The procedure begins with a Mutual Information affine registration of the image with the tissue probability maps. MI is computed from a 4x256 joint probability histogram.

See D'Agostino, Maes, Vandermeulen & P. Suetens. “Non-rigid Atlas-to-Image Registration by Minimization of Class Conditional Image Entropy”. Proc. MICCAI 2004. LNCS 3216, 2004. Pages 745-753.

Background voxels excluded Joint Probability Histogram

Background Voxels are Excluded

An intensity threshold is found by fitting image intensities to a mixture of two Gaussians. This threshold is used to exclude most of the voxels containing only air.

Spatially normalised BrainWeb phantoms (T1, T2 and PD) Tissue probability maps of GM and WM Cocosco, Kollokian, Kwan & Evans. “BrainWeb: Online Interface to a 3D MRI Simulated Brain Database”. NeuroImage 5(4):S425 (1997)

*

Image Registration

Figure out how to warp one image to match another * Normally, all subjects’ scans are matched with a common template

*

Current SPM approach

Only about 1000 parameters.

* Unable model detailed deformations

* *

A one-to-one mapping

Many models simply add a smooth displacement to an identity transform * One-to-one mapping not enforced Inverses approximately obtained by subtracting the displacement * Not a real inverse

Large deformation approximation Small deformation approximation

DARTEL

*

Parameterising the deformation

*

φ

(0)

(x) = x

1 *

φ

(1)

(x) =

∫

u

(

φ

(t)

(x)

)

dt

t=0 *

u

is a flow field to be estimated

Flow Field

Template

Initial Average Iteratively generated from 471 subjects Began with rigidly aligned tissue probability maps After a few iterations Used an inverse consistent formulation Final template

Grey matter average of 452 subjects – affine

Grey matter average of 471 subjects

Subject 1 Subject 2

VBM results of Patients with Alzheimer’s Disease

Both sets were smoothed with only 3mm FWHM Gaussian kernel, (FWE p>0.05)

Conventional SPM5 normalized image DARTEL normalized image

VBM results of Patients with Huntington’s Disease

Both sets were smoothed with only 6mm FWHM Gaussian kernel, (FWE p>0.05) 5 4 3 2 1 0

Conventional SPM5 normalized image

6 3 2 5 4 1 0

DARTEL normalized image

* * * * * * * *

References

Friston et al.

Spatial registration and normalisation of images.

Human Brain Mapping 3:165-189 (1995).

Collignon et al.

Automated multi-modality image registration based on

information theory. IPMI’95 pp 263-274 (1995).

Ashburner et al.

Incorporating prior knowledge into image registration.

NeuroImage 6:344-352 (1997).

Ashburner & Friston.

Nonlinear spatial normalisation using basis functions.

Human Brain Mapping 7:254-266 (1999).

Thévenaz et al.

Interpolation revisited.

IEEE Trans. Med. Imaging 19:739-758 (2000).

Modeling geometric deformations in EPI time series.

Neuroimage 13:903-919 (2001).

Ashburner & Friston.

NeuroImage in press (2005).

Ashburner . A fast diffeomorphic image reigstration algorithm. NeuroImage in press (2007).

Unified Segmentation.