Evaluation of 3D MRI Image Registration Methods

(1)

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2017

Evaluation of 3D MRI image

registration methods

(2)

Magnus Ivarsson LiTH-ISY-EX--17/5037--SE Supervisor: Andreas Robinson

cvl, Linköping University

Thobias Romu

amra

Magnus Borga

amra

Examiner: Maria Magnusson

cvl_{, Linköping University}

Computer Vision Laboratory Department of Electrical Engineering

(3)

Abstract

Image registration is the process of geometrically deforming a template image into a reference image. This technique is important and widely used within the field of medical IT. The purpose could be to detect image variations, pathological development or in the company AMRA’s case, to quantify fat tissue in various parts of the human body.

From an MRI (Magnetic Resonance Imaging) scan, a water and fat tissue image is obtained. Currently, AMRA is using the Morphon algorithm to register and seg-ment the water image in order to quantify fat and muscle tissue. During the first part of this master thesis, two alternative registration methods were evaluated. The first algorithm was Free Form Deformation which is a non-linear parametric based method. The second algorithm was a non-parametric optical flow based method known as the Demon algorithm. During the second part of the thesis, the Demon algorithm was used to evaluate the effect of using the fat images for registrations.

(4)

(5)

Acknowledgments

I would like to thank my supervisor Andreas Robinson and examinor Maria Mag-nusson at CVL for your great support during my master thesis. I would also like to thank my supervisors Thobias Romu and Magnus Borga for giving me the opportunity to complete my master thesis at AMRA and for your continuous support throughout the project.

Linköping, Januari 2017 Magnus Ivarsson

(6)

(7)

Notation

Notation Meaning IT Template image

IR Reference image

Ω Domain where an image I is defined u Displacement field

δu Incremental displacement field c Certainty field

∇_I _{Image gradient}

D Distance/Similarity measure

T Transformation (unless something else is specified) R Regularization

x Bold lowercase character indicate vector A Bold uppercase character indicate matrix h_{I, Ii} _{Scalar product between two images} ∇_{· u} _{Divergence of a displacement field}

∗ _Convolution

T ◦ I Transformation applied to an image

Abbreviation Meaning

MRI Magnetic resonance imaging VAT Visceral adipose tissue

ASAT Abdominal subcutaneous adipose tissue SAT Subcutaneous adipose tissue

LULB Left upper leg back (muscle) LULF Left upper leg front (muscle) RULB Right upper leg back (muscle) RULF Right upper leg front (muscle)

IP In phase (water + fat image)

(10)

(11)

1

Introduction

This chapter is intended to give the reader the necessary background information and motivation of why image registration of MR images is useful and the aim of the master thesis. The problem formulation and the limitations of this thesis are also described.

1.1 Background

With an aging population, an efficient and improved health care is essential to handle the demands we put on our health care system. Medical imaging and image registration of medical images will play an important role in the future. In today’s society, health and metabolic status for individuals are often obtained by indirect measurements such as BMI and waist circumference. However, alternate methods exist: AMRA is a company that specializes in precise body composition measurements that can improve and customize treatments for individuals with high metabolic risk. This is done by generating a water and fat image from an MRI scan followed by a quantification of fat and muscle tissue in different regions of the body. This quantification is possible because of two reasons: First, the water image indicate which voxels that contain muscle tissue while the fat image indicate which voxels contain fat tissue. Second, the images are normalized so that the voxel intensity correspond to the concentration of fat or muscle tissue within that particular voxel. An example of a water and fat image can be seen in figure 1.1 and images of segmented adipose and muscle tissue can be seen in figure 1.2.

(12)

Figure 1.1: The image to the left illustrates a water image, the image in the middle illustrates a fat image and the image to the right illustrates an IP image (water + fat) in the coronal plane. Note that the water signal is strong where the fat signal is weak and vice versa. Image source: AMRA.

Figure 1.2: The images to the left illustrate the upper leg muscle regions. Blue pixels represent left upper leg back, yellow represents left upper leg front, magenta represents right upper leg back and cyan represents right upper leg front. In the image to the right, red pixels represent VAT and blue pixels represent ASAT. Image source: AMRA.

(13)

1.1 Background 3

Examples of measurements AMRA provide are

• VAT - Visceral Adipose Tissue (intra-abdominal fat)

• ASAT - Abdominal Subcutaneous Adipose Tissue (pincheable fat) • Thigh Muscle Volume

• Lean Muscle Tissue Volume • Total Adipose Tissue

• Additional Muscle Group Volumes

The reason why the body is segmented into different regions is that it often mat-ters where fat tissue is located. For instance, research has shown that large vol-umes of visceral adipose tissue is connected to diabetes, liver disease, and cancer [15]. On the other hand, subcutaneous fat is much less likely to have adverse health effects [16]. Tools that are used today, like BMI and waist circumference, are well suited to determine health and metabolic status on a population basis but can sometimes be crude when used to determine the body composition of an individual. An image that illustrates the difference between BMI and VAT can be seen in figure 1.3.

Figure 1.3:Example of six men with a BMI of 21 where their body composi-tion and metabolic risk have a high variacomposi-tion. The blue color represents ab-dominal subcutaneous adipose tissue and red color represents visceral adi-pose tissue. Image source: AMRA.

(14)

So far, the reasons why precise body composition measurements are necessary have been explained. However, the method by which the regions of interest are found, i.e. how to segment the image, has not yet been described. A simple but time-consuming way to segment an image is manual labelling of each voxel. An-other approach, which is more cost effective, is to use image registration methods to automatically segment an image with the aid of a computer. Note that the precision is not guaranteed to be better by manual segmentation compared to au-tomatic. Today, AMRA is using the Morphon algorithm in order to register and segment an MR image.

1.1.1 Aim

The aim of the master thesis is to evaluate different registration algorithms and to determine whether the results given by the Morphon can be improved. The evaluation was carried out as explained in section 2.5.1 and 3.6.

1.2 Problem Formulation

The master thesis has been divided into two parts. The first part was devoted to answer the following questions:

• Is it possible to obtain better segmentation with the parametric method Free Form Deformation (explained in section 2.1.1) compared to the Morphon? • Is it possible to obtain better segmentation with the Demon algorithm

(ex-plained in section 2.4.2) compared to the Morphon?

In the second part, the Demon algorithm was used to answer the following ques-tions:

• Is it possible to obtain better segmentation if the fat image is used as input to the registration algorithm?

• Is it possible to obtain better segmentation if a combination of the fat and water image is used as input to the registration algorithm?

1.3 Limitations

Because of the limited time frame of the master thesis, the methods chosen for evaluation are within the group of parametric and non-parametric methods. Due to the high resolution of the 3D images and the computational demands of the algorithms, only a modest number of images were used for evaluation. How-ever, the number of samples is still enough to indicate whether an algorithm has a lower, higher or the same potential than the Morphon algorithm.

(15)

1.3 Limitations 5

Since the focus of this thesis is image registration, the theory behind MR imaging and how the water and fat images is obtained are not within the scope of this thesis. Neither is the theory behind the normalization of the images to do a quan-titative measurement of fat and muscle tissue. The interested reader is referred to [11], [10] and [17].

(16)

(17)

2

Theory

This chapter primarily explains the theory behind image registration but also describes some metrics that can be used to evaluate a registration.

2.1 Image registration - an overview

This section provides an overview of the most common image registration tech-niques used today. The content given in this section is based on [3].

Before presenting different types and groups of algorithms, it is important to un-derstand the main concept. The idea is to find a transformation T that deforms a template image IT so that the deformed template image T ◦ IT is similar to a

reference image IR. The transformation T can be a transform applied to all voxels

in the image IT, which is referred to as a parametric transform. It can also be a

non-parametric transform where all voxels are allowed to be displaced separately. In the latter case, the transformation is seen as a displacement field which is de-noted by u instead of T . These types of transformations are further explained in section 2.1.1 and 2.1.2. An example of a template image before and after it is deformed into a reference image can be seen in figure 2.1.

All registration algorithms, no matter which group they belong to, aim to mini-mize some kind of cost function described by

(T ) = D(IR, IT, T ) + R(T ), (2.1)

where D is a distance measure and R is a regularization term. The term D is used to determine how different the reference image IRand the template image

IT are after the transformation T has been applied to IT. The term R is added

(18)

to regularize the solution so that the displacement field becomes homogeneous and to penalize unreasonable transformations. Some methods, however, are not constructed as a pure optimization problem. Instead the problem is divided into two parts:

1. Finding the optimal transform T by a distance measure. 2. Regularizing the transform.

The main reason for dividing the problem into two parts is to simplify it and facil-itate the parameter tuning. Non-parametric image registration methods, which are explained in section 2.1.2, typically do this separation.

Figure 2.1:This figure illustrates the effects of an image registration with the water images in the axial plane of a body. The image to the left is a template image IT, the image in the middle is a deformed template image T ◦ IT and

the image to the right is a reference image IR. Note that the template image is

more similar to the reference image after the deformation than before. Image source: AMRA.

2.1.1 Parametric methods

In parametric methods, the aim is to find a set of parameters p that minimizes a cost function under the constraint of a certain type of transformation. The idea of using a parametric instead of a non-parametric transform is that the parametriza-tion explicitly regularizes the soluparametriza-tion [13]. First, linear transformaparametriza-tions will be explained based on [3]. Second, non-linear transformation is briefly covered based on [13].

Linear transformation

Examples of linear transformation types are • Rigid - translation and rotation. • Similarity - rigid and scaling. • Affine - similarity and skewing.

(19)

2.1 Image registration - an overview 9

The general form of these transformations can be described as

Tp(x) = x +

X

k

pkBk(x), (2.2)

where pk is a scalar, x is the coordinate of a voxel, Bk is a multivariate basis

func-tion and Tp(x) is the transformed coordinate. The generic cost function described

by (2.1) can then be modified into

(p) = D(IR, IT, T (p)) + R(p). (2.3)

In general, the regularizing term R is omitted since the idea behind a parametric transform is that the constraints of the parameters themselves will regularize the solution.

A rigid transform in 2D for instance, which consists of a rotation R and a transla-tion t, can be written as

Tp(x) = Rx + t ="p1 −_p₂ p2 p1 # "x y # +"p3 p4 # ="x −y 1 0 y x 0 1 #             p1 p2 p3 p4             = B(x)p (2.4)

where the columns of B(x) are the basis functions Bk(x).

Non-linear transformation

In order to allow for more complex deformations, non-linear transformations are the natural step after the affine transformations. Note that non-linear transforma-tions should not be confused with non-parametric methods which are explained in section 2.1.2. Non-linear transformations are still governed by a set of pa-rameters while non-parametric transformations allow each voxel to be displaced separately.

Thin-plate spline (TPS) is a well known non-linear and grid-based registration algorithm. The idea is to place an evenly spaced grid of control points in the template image and find an affine transformation for these points. The transfor-mation for the remaining voxels is then given by,

T (x) = Ax + t +

n

X

i=1

(20)

where A is an affine transformation matrix, t a translation vector, n the number of control points, xi is the coordinate of the ithcontrol point and wiis the weight

of the ithcontrol point. In TPS, all control points are included in each

transforma-tion i. e. they have global support. This makes the transformatransforma-tion more robust but do not allow as complex local deformations as Free Form Deformation which will be explained in section 2.4.3.

2.1.2 Non-parametric methods

As described in section 2.1.1, parametric methods aim is to find a transformation T that deforms the template image given certain constraints. For non-parametric methods, the deformation is no longer bound to a certain type of transformation. Instead, all voxels of the template image are displaced separately which allows for very complex deformations. Since the deformation is no longer a product of a transformation T applied to all voxels, the transformation is replaced by a displacement field u. This displacement field u contains the displacement of all voxels within the image.

Two ways of finding the displacement field for a non-parametric method will be covered in section 2.4 where the Morphon and Demon algorithms will be ex-plained. These are two well known non-parametric image registration algorithms where Morphon make use of the local phase to find the displacement while De-mon is based on optical flow.

The content in this section is based on [3] unless anything else is specified. Field Accumulation

To be able to find a large and accurate deformation, an iterative approach through scale-space where the displacement field is accumulated is often preferred. There are many ways to accumulate the displacement field and this thesis explains three alternatives, see (2.6), (2.7) and (2.8).

The simplest way is to add the incremental displacement δu to the current accu-mulated displacement u by

u ← u+ δu. (2.6)

However this approach can cause unwanted tears and folding which is not phys-ically possible for a deformed body. The usage of composite accumulation is preferred in this situation,

u ← u ◦δu = u(x + δu(x)) + δu(x). (2.7) Another way to avoid tears and foldings of the displacement field is by accumu-lating u as

(21)

2.1 Image registration - an overview 11

u ← u ◦exp (δu), (2.8)

which ensures that the displacement field stays diffeomorphic. This means that the field is invertible and that both the field and its inverse are smooth or, to put differently, that its Jacobian determinant is positive.

Field Regularization

An important concept in image registration is the displacement field regulariza-tion. In section 2.3.1 the regularization of the displacement field was a part of the error function explicitly. However in non-parametric image registration methods the incremental and/or the accumulated displacement fields need to be regular-ized to make the problem well-posed. The simplest and most straight-forward way is to convolve the field with an isotropic low pass filter

u ← u ∗giso. (2.9)

This filter will efficiently smooth the displacement field in homogeneous regions, but it will also smear out edges and lines which is not preferred. This can be mitigated by adaptive filtering where the filter size is modified depending on the local region.

The following content describes a method to perform anisotropic filtering sug-gested by [4]. The examples and equations are given in 2D but can be extended to 3D. The first step is to compute the local structure tensor T of the displacement field and do an eigenvalue decomposition

T =

2

X

k=1

λkˆekˆeTk. (2.10)

A large eigenvalue λk indicates that there is a distinct structure orthogonal to

their corresponding eigenvector ˆek. Therefore low pass filtering should only be

executed in the directions of the eigenvectors that correspond to small eigenval-ues. There are basically three different cases where:

1. No apparent structure exists, both λ1and λ2small.

2. Distinct structure exists in one direction, λ1large and λ2small.

3. Distinct structure exists in two orthogonal directions, λ1and λ2large.

(22)

C =

3

X

k=1

γkˆekˆeTk, (2.11)

where γk is mapped from λk such that 0 ≤ γk ≤ 1. Further, an adaptive filter

gadapcan be created by a linear combination of three filters:

• glarge - isotropic filter which suits case 1.

• ganiso- anisotropic filter where the orientation is governed by ˆe1which suits

case 2.

• gsmall- small isotropic filter which suits case 3.

The adaptive filter is given by

gadap= ||C|| · ((1 − γ2 γ1 ) · ganiso+ γ2 γ1 · gsmall) + (1 − ||C||) · glarge. (2.12)

The parameters γkcan be found by

T1p= T ∗ g, (2.13) γ1 = m(||T1p||, σ , α, β), (2.14) γk = γ1 k Y l=2 µ( λl λl−1 , αl, βl), k = 2, ..., N , (2.15)

where g is a Gaussian low pass filter and N is the number of dimensions. The m and µ functions are used to control the transition between case 1-2 and between case 2-3. Examples of the functions m and µ and the purpose of the parameters σ , α and β are not given by [4].

2.1.3 Landmark based registration

Landmark based image registration is a sparse method that is based on anatomi-cal landmarks, which are a set of interest points that can be found in most human beings. An image that illustrates a set of landmarks and their correspondence between two images can be seen in figure 2.2. The content in this section is pri-marily based on [12] unless anything else is specified.

The idea is to find a transformation T that minimizes the distance between the landmarks in the template and reference image. Since the correspondence be-tween one point in the template image and one point in the reference image is known, there is no need of introducing a similarity measure. Which transfor-mation that minimizes the distance between the landmarks depends on the con-straints of the transformation. As explained in section 2.1.1, such concon-straints

(23)

2.2 Distance measures 13 could be to allow a translation of the points or a rigid, similarity, affine or non-linear transformation.

An important detail when describing landmark based registration is how to actu-ally find the coordinates of the landmarks. The simplest method is to manuactu-ally locate each landmark in the images. This can be an accurate but time consuming approach that is not easily scalable. To overcome the problem of manual place-ment, a classifier can be trained to find each landmark which would make this approach scalable. An example of how a convolutional neural network can be trained to detect landmarks on distal femur in 3D MR images is presented in [23].

Figure 2.2:Example of correspondences between landmarks in two images. The green circles indicate the landmarks and the red lines in-between them determine their correspondence.

2.2 Distance measures

To determine whether a transformation applied to IT makes it similar to IR

de-pends primarily on which distance measure is used. This section gives an intro-duction to the most common distance measures used in image registration based on [12], [3] and [13] unless anything else is specified. Some measures are simple, have a low computational cost, while other methods can handle more complex contexts but comes at a higher computational cost. The transformation in this section will be regarded as a displacement field, u.

Before going through the different methods, the topic of multi-modality has to be mentioned since some distance measures can handle this situation better than others. Multi-modality refers to images generated by different settings of a MRI

(24)

scan or different scanning methods (ex MRI and CT), which vary in intensity and tissue contrast. To name one example, it can be an additive or a multiplicative intensity offset between the images. This makes image registration a much more complex problem since simple methods, which assume that similar structures have similar intensity values, won’t work in this case. In figure 2.3 a T1 and T2 weighted MRI scan can be seen. Note how they have similar shapes and struc-tures, but completely different intensity levels in some areas. Since the focus of this master thesis is on image registration, T1 and T2 weighted MRI scanning will not be explained. It should also be mentioned that this thesis only considered mono-modal image registration.

Figure 2.3:The image to the left is a T1 weighted MR image and the image to the right is a T2 weighted MR image. Image source: Wikipedia [7].

2.2.1 SSD

One well known intensity based method is sum of squared differences (SSD), DSSD(IR, IT, u) = Z Ω IR(x) − IT(x + u(x)) 2 dx, (2.16)

where Ω is the image domain, x is a voxel coordinate and u(x) is the displacement for that particular voxel.

This distance measure assumes that the intensities match when the correct trans-formation is applied to the template image. This method work well for mono-modal image registration but not for multi-mono-modal registration since it is intensity based.

(25)

2.2 Distance measures 15

2.2.2 Normalized Cross Correlation

The cross-correlation CC between a template image IT and reference image IRis

defined as

CC(IR, IT, u) = < IR, IT(u) > =

Z

Ω

IR(x)IT(x + u(x))dx. (2.17)

Normalized cross-correlation is given by

N CC(IR, IT, u) = < IR, IT(u) >

< IR, IR> < IT(u), IT(u) >

, (2.18)

and the corresponding distance measure is defined as

DN CC(IR, IT, u) = 1 − N CC(IR, IT, u)2. (2.19)

This approach can handle images where the intensities are related by a multi-plicative factor because of the normalization in (2.18).

2.2.3 Normalized Gradient Field

The normalized gradient field, N GF, was first presented in [6] and the following content is based on that article. The equation for N GF is given by

N GF(I) = √ ∇I

∇_IT∇_{I +}2, (2.20)

where ∇I is the image gradient and is a parameter which makes the distance measure robust against noise. It is included to indicate the smallest change in the image I that should be interpreted as a distinct edge. This parameter is given by

 = η V

Z

Ω

|∇_I(x)|dx _(2.21)

where η is the estimated noise level in the image and V is the volume of the image domain. The corresponding distance measure to NGF is given by

DN GF(IR, IT, u) =

Z

Ω

1 −N GF(IR(x))TN GF(IT(x + u(x))

2

dx. (2.22)

There are two reasons why this is a powerful distance measure. First, the gradient image is normalized which solves the problem with a multiplicative intensity difference between the images. Second, since gradients are computed by voxel differences, any additive intensity difference will be canceled by that operation.

(26)

2.2.4 Mutual Information

Mutual information, MI, is a distance measure based on statistical information rather than spatial information [5]. The distance measure is given by

DMI(IR, IT, T ) = Z R pRlog (pR)dr + Z R pT log (pT)dt − Z R2 pR,T log (pR,T)drdt, (2.23) where pR and pT are the marginal densities and pR,T is the joint density over

intensity values. This distance measure is the only one which can handle a non-linear mapping of intensities presented in this thesis.

2.3 Image registration as an optimization problem

As mentioned in section 2.1, image registration aims to minimize a cost function under certain constraints. How this cost function and constraints are formulated may vary between applications and will heavily affect the end result. A naive approach would be to let all voxels be displaced separately and to use SSD from (2.16) as a distance measure

(u) = DSSD(IR, IT, u) =

Z

Ω

(IR(x) − IT(x + u(x)))2dx. (2.24)

In this case, the error function would be easy to minimize and the cost would probably be low. However, the displacement field would not be particularly ac-curate or unique.

2.3.1 Regularization methods

In order to improve the probability of finding an accurate and plausible displace-ment field, i.e. making the problem well-posed, a regularization term R can be added to the cost function

(u) = DSSD(IR, IT, u) + R(u). (2.25)

This section will briefly explain the most common regularization terms, i.e. dif-fusion, curvature, elastic and fluid based on [3], [13], [2] and [12].

Diffusion

Diffusion regularization aims to penalize the inhomogeneities in the displace-ment field by, in each dimension, taking the sum of the norm of the gradient of u. The equation is given by

(27)

2.4 Methods 17 R(u) = 1 2 Z Ω d X i=1 ||∇_u_i||2_dx _(2.26)

where i indicate the dimension of u and d is the number of dimensions. Note that (2.26) is equal to zero if the transformation corresponds to a translation.

Curvature

Contrary to diffusion regularization, curvature regularization allow more com-plex deformations by taking the Laplacian instead of the gradient of the displace-ment field. As long as the deformation is affine, the regularization term will be equal to zero. The equation is given by

R(u) = 1 2 Z Ω d X i=1 (∆ui)2dx, (2.27)

where ∆ is the Laplace operator. Elastic

Elastic regularization has a physical interpretation where an elastic material’s energy potential is minimized. The equation is given by

R(u) = 1 2

Z

Ω

µ < ∇u, ∇u > + (λ + µ)(∇ · u)2dx, (2.28)

where λ and µ are the so-called Lamé constants. Fluid

The equation for fluid regularization is the same as for elastic regularization ex-cept that it is applied to the velocity field, e.g. the incremental displacement field δu.

2.4 Methods

This section describes the general theory behind the methods that were evaluated in this thesis. The first two methods, Morphon and Demon, are both well known non-parametric methods. The last method, Free form deformation, is a typical non-linear and grid based parametric method.

(28)

2.4.1 Morphon

The Morphon algorithm was first introduced [9] and the content in this section is based on that article and [3]. One of the ideas is to estimate the local displace-ment u by computing the phase difference between the template and reference image. This can be done by convolving both images with a set of N quadrature filters {fk}Nk=1in N different directions { ˆnk}Nk=1. Consequently the filter responses,

{_q_T

k}

N

k=1and {qRk}

N

k=1, for the template and reference image are given by

qTk = IT ∗fk, (2.29)

qRk = IR∗fk, (2.30)

where fk is a quadrature filter. In practice, fk is represented by two filters, one

filter that represents the real part and one filter that represents the imaginary part of the complex filter. An illustration of how the phase of a quadrature response is related to different types of edges can be seen in figure 2.4.

Figure 2.4: Illustration of the phase ϕ from a quadrature response. Image source: [1]

The conjugate product of the filter responses can be denoted as Qk = qTkq

∗

Rk, (2.31)

and the phase differences pkcan then be computed by taking the argument of the

conjugate product of the filter responses

pk = arg(Qk). (2.32)

In order to increase the robustness against noise, the local phase differences can be weighted by the magnitude of the conjugate product

(29)

2.4 Methods 19

ck =

p

|_Q_k|_, _(2.33)

where c stands for certainty. A clear edge, for instance, will result in a large Qk

while noise will give a smaller response. The certainties are used to weight the displacements u when estimating the displacement field (2.35). They can also be further regulated by incorporating the phase difference,

ck= |Qk|1/2cos2(pk/2), (2.34)

where small differences in phase will give a higher certainty. The simplest way to estimate the displacement, which was presented in [19], is by a weighted sum given by u= PN k=1ckpknˆk PN k=1ck , (2.35)

where ck is the certainty and pk is the phase difference in the direction of ˆnk.

As explained in section 2.1.2 and suggested in [9], the displacement field should be accumulated incrementally while iterating from coarse to fine scale. The accu-mulated displacement u and certainties c are then updated with

u= cu+ δc(u + δu)

c+ δc , (2.36)

c=c

2_{+ δc}2

c+ δc , (2.37)

where δc is the incremental certainty and δu is the incremental displacement on a specific scale. To update the displacement field with help of the certainty map characterizes the Morphon algorithm and is one of the benefits of using the local phase as a distance measure. Another advantage of the local phase difference is its invariance to any intensity differences between the template and reference im-ages.

Since the registration problem is ill-posed and there exists multiple displace-ments, the field can be regularized by normalized convolution where the certain-ties are incorporated,

δu = (δuδc) ∗ g

δc ∗ g , (2.38)

(30)

2.4.2 Demon

The Demon algorithm, first described in [21] and [20], was inspired by the con-cept of optical flow to find a displacement u of the template image. The content in this section is based on those articles and [3].

Optical flow refers to the motion between two images taken with a small time difference ∆t apart. Assuming that the brightness in the two images is constant, the following equation must hold,

I(x, t) = I(x + ∆x, t + ∆t). (2.39) A first order Taylor expansion of (2.39) yields

∇_IT_∆x_{= I}_t_, _(2.40)

where It is the derivative of I with respect to t and ∇I is the spatial gradient of

the image. Since the solution to (2.40) is not unique, there exists many solutions where the simplest solution is given by

∆x= It

||∇_I||2∇I, (2.41)

which also served as an inspiration to the Demon algorithm. Since Itdenotes the

change of time and IR−IT can be seen as two images taken one time step apart,

Itcould be replaced by IR−IT and ∆x by u in (2.41) which gives

u= IR−IT

||∇_I_R||2∇IR. (2.42)

A modified version, which is more stable for small values of ∇IRis given by

u= IR−IT ||∇_I_R||2_{+ (I}

R−IT)2

∇_I_R_. _(2.43)

Ref. [20] suggests the use of a multiscale implementation, similar to the Morphon algorithm, for several reasons. First of all it speeds up the implementation, the additional computational cost of performing a computation on all scales is 1/3 of the computational cost for the finest scale. In addition, the convergence is faster. Finally, macroscopic features are more stable for human anatomy.

Since (2.43) is not a unique solution, it is necessary to regularize the incremental displacement field δu with a Gaussian low pass filter,

δu = δu ∗ g, (2.44)

(31)

2.4 Methods 21

2.4.3 Free Form Deformation

Free form deformation - FFD - is a well known non-linear and grid based method similar to Thin-plate splines described in section 2.1.1. The following content is based on [18] and [13].

The name Free Form Deformation is slightly misleading since the deformation is still governed by a parametric transformation under a few constraints, as op-posed to non-parametric deformations. The idea is to place an evenly spaced mesh (nx×ny×x × nz) of control points in the template image and find an affine

transformation for these points. The deformations in between these points are found by non-linear interpolation using 1D cubic B-splines. The deformation of a single voxel is given by

T (x, y, z; p) = 3 X k=0 3 X m=0 3 X n=0 Bk(ux)Bm(vy)Bn(wz)pi+k,j+m,l+n, (2.45)

where x, y and z are the coordinates of a voxel and i = bx/nxc−1, j = by/nyc−1, l =

b_z/n_zc −_{1 and p are the control points. The local B-spline coordinates (u}_x_{, v}_y_{, w}_z₎ are given by

ux= x/nx− bx/nxc, (2.46)

vy= y/ny− by/nyc, (2.47)

wz= z/nz− bz/nzc, (2.48)

and the B-spline functions are given by

B0(u) = (1 − u)3/6, (2.49)

B1(u) = (3u3−6u2+ 4)/6, (2.50)

B2(u) = (−3u3+ 3u2+ 3u + 1)/6, (2.51)

B3(u) = u3/6. (2.52)

There are some similarities but also some differences between FFD and TPS. Both methods aim to find an affine transformation of the control points. The difference, however, is by which method the displacements between the control points are interpolated and how many of the control points affect each interpolation. In TPS, the deformation of all control points impact the deformation of all other voxels which means that they have global support. In FFD, the deformation of a voxel is only affected by the nearest control points, which means that they have local sup-port. The consequence is that FFD allows for more complex local deformations while TPS tends to be more robust.

(32)

2.5 Evaluation

In order to determine how well an image registration method works, some kind of similarity metric has to be used. The first type of metric, body composition measurements, evaluates the effect of the image registration. The second type, segmentation metrics, evaluate how well a deformed mask pertaining to a tem-plate image overlap with the ground truth mask pertaining to the reference im-age.

2.5.1 Body Composition Measurements

Given a registered and segmented image, it is possible to compute the body com-position measurements mentioned in section 1.1. The muscle measurements be-low are given by [8]. The muscle volume is given by

MV = X

∀_voxels

Mmuscle· Mbody· Vvoxel, (2.53)

where Mmuscleis a deformed binary muscle mask pertaining to the template

im-age, Mbody a binary body mask, Vvoxel the volume of a voxel connected to the

reference image. The fat free muscle volume MVFFis given by

MVFF=

X

∀_voxels

(f < 0.5) · Mmuscle· Mbody· Vvoxel, (2.54)

where f is the fat-concentration in a particular voxel connected to the reference image. This means that each voxel within Mmuscle with a fat concentration less

than 0.5 is considered to be a fat-free muscle voxel.

Estimation of VAT and ASAT can be done in a similar manner,

V AT = X

∀_voxels

f · MV AT· Mbody· Vvoxel, (2.55)

ASAT = X

∀_voxels

f · MASAT· Mbody· Vvoxel, (2.56)

where MV AT and MASAT are the binary VAT and ASAT masks.

Bland-Altman plot

As a complement to the standard statistics like the mean µ, standard deviation σ , median xmed, minimum xmin, maximum xmaxand range xrangeof a sample set,

a Bland-Altman plot can be used to determine the relation between two differ-ent methods. Consider two sets of n samples (e.g. n differdiffer-ent image volumes), {_s1

(33)

2.5 Evaluation 23

method and the second set contain the measurements of the second method. For the ithmeasurement, the average siav and difference sdif fi between the two

meth-ods are given by

s_avi = s i 1+ s2i 2 , (2.57) si_{dif f} = s₁i −_si 2. (2.58)

The sample mean difference also known as the bias ˆµdif f and the sample

stan-dard deviation of the difference ˆσdif f are given by

ˆ µdif f = 1 n n X i=1 si_{dif f}, (2.59) ˆ σdif f = v t 1 n − 1 n X i=1 s_{dif f}i −_µ_ˆ_{dif f}_. _(2.60)

In the Bland-Altman plot, the difference {s1

dif f, ..., s n

dif f}is plotted against the

av-erage {s1_av, ..., sn

av}as a scatter plot for the n measurements. In addition, the mean

difference ˆµdif f is visualized as a straight line and the standard deviation of the

difference ˆσdif f is used to visualize the 95% confidence interval assuming that

the difference follows a Gaussian distribution. An example of a Bland-Altman plot can be seen in figure 2.5.

(34)

Figure 2.5:The green dots represent each sample in a dataset. The samples are located depending on the average (2.57) and difference (2.58) between the two methods. The blue line (bias) is given by (2.59). The red dotted lines are given by ˆµdif f ±1.96 ˆσdif f which represent the 95% confidence interval

for a Gaussian distribution.

2.5.2 Segmentation metrics

To determine how accurately the algorithms register a reference image, a few well known similarity metrics are presented in this section. These metrics are not, by definition, specific to image registration and segmentation, but are applicable in this context. To understand the similarity metrics, the concept of true positives (TP), false positives (FP), true negatives (TN) and false negatives (FN) have to be introduced. In the field of 3D image segmentation they can be defined in the following ways:

• True positives, TP, refer to the sum of voxels which are correctly labelled as foreground.

• False positives, FP, refer to the sum of voxels which are falsely labelled as foreground.

• True negatives, TN, refer to the sum of voxels which are correctly labelled as background.

• False negatives, FN, refer to the sum of voxels which are falsely labelled as background.

(35)

2.5 Evaluation 25

The number of true positives or false negatives in an image provide, by them-selves, very little information. However the relationship between these measure-ments do. To name an example, if the number of true positives in an image is 500, is that a good or bad result? If the number of false positives is 4 and the number of false negatives is 6, it is reasonable to say that the answer to that question is yes. On the other hand, if the number of false positives is 400 and the number of false negatives is 600, it is likely that the answer is no.

There are many ways to relate these labels to each other and the following content introduces a few ways to determine similarity between a deformed template im-age and a reference imim-age. In order to fully understand the concept of TP, FP, TN, FN and the similarity measurements described below, an image that illustrates these concepts is provided in figure 2.6.

Precision

Precision is a measure which describes how accurately a method determines true positives within the set of all attempts to do so. When a method gives a positive answer, the precision measurement informs us of how likely the answer is to be correct. The equation is given by

P recision = T P

T P + FP. (2.61)

Images that illustrate how the precision measurement is affected in different situ-ations can be seen in figure 2.6.

Recall

Recall is a measure which describes how likely a method is to determine all true positives, regardless of false positives. The equation is given by

Recall = T P

T P + FN. (2.62)

Images that illustrate how the recall measurement is affected in different situa-tions can be seen in figure 2.6.

(36)

Figure 2.6: This image illustrate the back and front leg muscles in the coronal plane. The yellow pixels represent true positives, magenta repre-sent false positives and cyan reprerepre-sent false negatives. In the left image, P recision < 1.0 and Recall = 1.0. In the right image, P recision = 1.0 and Recall < 1.0. Image source: AMRA.

Specificity

Specificity is a measure which describes how accurate a method determines true negatives within the set of all attempts to do so. When a method gives a true negative answer, the specificity measurement informs us of how likely the answer is to be correct. The equation is given by

Specif icity = T N

T N + FP. (2.63)

Dice

A common similarity metrics explained in [22] is the Dice Similarity Coefficient given by DSC = 2T P 2T P + FP + FN = 2 · ||ΩI_R∩ ΩI_T|| ||Ω_I_R||_{+ ||Ω}_I T|| , (2.64)

where ΩIR and ΩIT are the regions of the template and reference image. The

advantage of Dice compared to precision and recall is that both false positives and false negatives are included in the metric. In the images illustrated in figure 2.6, the Dice Similarity Coefficent will accurately describe the relation between true positives and the sum of false positives and false negatives. The downside

(37)

2.5 Evaluation 27

of Dice is that it doesn’t describe which type of error, false positives or false neg-atives that occurs the most.

(38)

(39)

3

Method

This chapter describes the methods used to answer the problems formulated in section 1.2. First, the Free Form Deformation (FFD) and Demon algorithms were compared to the Morphon algorithm. Second, the Demon algorithm was further investigated.

Section 3.1 describes the exact variants chosen for the Morphon, FFD and Demon algorithms.

Section 3.2 describes the further investigations with the Demon algorithm. Section 3.3 describes the regions that are of interest to segment.

Section 3.4 describes the datasets and the properties of the templates (prototypes) and reference images.

Section 3.5 describes how the prototypes were registered to the reference image, how the most appropriate prototypes are selected and how the segmentation of different tissue are performed.

Section 3.6 explains the evaluation and its corresponding measurements.

3.1 Algorithm comparison

Before presenting the algorithms, it is important to mention that FFD and De-mon both measure voxel intensity differences. It was possible to use this distance measure because the images were normalized as mentioned in section 1.1. All algorithms evaluated in this thesis used a multiscale approach where the fi-nal displacement was found by iterating the algorithms several times on several scales and accumulating the displacement field accordingly. When describing an algorithm, a table is included which contains the parameters for the particu-lar implementation used in this thesis. The parameter downsampling factor is a

(40)

scalar which determines how much the images were downsampled for each scale. Start scale determines how many times the images were downsampled were the first iterations were made and finally the end scale determines on which scale the last iterations were made. The finest scale, which is the original size of the image, is equal to 0. In order to avoid any confusion, an example is provided:

Assume that the start scale is 3, the end scale is 0 and the downsampling factor between one scale and the next is 2. The scale were the first iterations are made are generated by downsampling the original image by 23= 8 and the second scale were generated by downsampling the original image by 22 = 4 and so on. The final scale were generated by downsampling the original image by 20 = 1, which is the finest scale and the original image itself.

A flowchart which describes the different steps in the registration process can be seen in figure 3.1. The results from this part are found in section 4.1.

(41)

3.1 Algorithm comparison 31

Figure 3.1:A flowchart that describes the different steps in the registration process.

(42)

3.1.1 Morphon

The Morphon was regarded as the baseline method of this thesis and was in-cluded to evaluate the results of FFD and Demon. The Morphon algorithm is implemented by AMRA in Python and constructed as described in section 2.4.1 where the certainties where computed as in (2.34). The specific parameters of their implementation are not specified in this thesis.

3.1.2 Free form deformation

Section 2.4.3 explains the general idea of FFD and how to find the transformation in-between the control points. This section describes the method used to find the displacements of the control points.

As mentioned in section 2.4.3, the idea of the FFD algorithm is to find an affine transformation of the control points and then interpolate of the displacement for the remaining voxels with cubic B-splines. In order to find an accurate transfor-mation which remain affine for the control points, the problem was formulated as a minimization problem with SSD (2.16) as a distance measure and curvature regularization (2.27) of the control points.

SSD was used since the images are normalised and any other distance measure would only cost more in terms of computational speed. As explained in section 2.3.1, curvature regularization implicitly adds a cost for non-affine tions of the control points. This does not necessarily mean that the transforma-tion will be affine, only that the regularizatransforma-tion term penalizes non-affine transfor-mations.

The cost function is given by

min u1,...,n(u1,...,n) = 1 2 Z Ω (IR(x) − IT(x + u(x))2dx + λ ·1 2· n X i=1 (∆ui)2, (3.1)

where u(x) is the displacement of a voxel, λ is a weight that controls the regu-larization, ui is the displacement of the ithcontrol point and n is the number of

control points.

The optimal solution to (3.1) was found with the implicit Euler method. A list of the parameters used for the FFD algorithm in this thesis are given by table 3.1. The implementation used for the FFD algorithm is a part of the Medical Image Registration Toolbox (MIRT) created by Andriy Myronenko in Matlab [14].

(43)

3.1 Algorithm comparison 33 Table 3.1:Parameters for the FFD algorithm

Parameter Value

Downsampling factor 2

Start scale 2

End scale 0

Iterations per scale [100 100 100] Control point distance 8 voxels

λ 0.01

3.1.3 Demon

The Demon algorithm that was evaluated had the following specifics. The incre-mental displacement field estimate δu was given by

δu = IR−IT ||∇_I_R||2_{+ (I}

R−IT)2

∇_I_R_, _(3.2)

explained in section 2.4.2. The gradient was computed by the MATLAB®function gradientwhich gives the numerical gradient. The displacement field u and the incremental displacement field δu were accumulated with a simple sum,

u ← u+ δu, (3.3)

explained in section 2.1.2. The displacement field u and the incremental displace-ment field δu were regularized with normalized convolution,

u= (uc) ∗ g

c ∗g , (3.4)

δu = (δuc) ∗ g

c ∗g , (3.5)

where g is a Gaussian low pass filter. The certainty field is given by the gradient magnitude of the reference image,

c= |∇IR|, (3.6)

which is constant for each scale. A list of the parameters used for the Demon algorithm are given by table 3.2. The implementation of the Demon algorithm is a part of the Fordanic/Image-registration Toolbox created by Daniel Forsberg in Matlab [3].

(44)

Table 3.2:Parameters used for the Demon algorithm. Note that the last value of the parameter iterations per level is the number of iterations for the end scale.

Parameter Value

Downsampling factor 2

Start scale 3

End scale 0

Iterations per scale [50 50 50 25] Gaussian filter size 5 x 5 x 5 Gaussian filter sigma 1 voxel

3.2 Further investigation

In the first part of the thesis, the water images were used to register the proto-types to the reference images for each algorithm. In this part, the effect of in-cluding the fat images in the registrations was evaluated. This modification is applicable to all the algorithms presented so far and the Demon algorithm (with the water images as input) was chosen as the basis algorithm since it is the fastest and easiest to modify.

The first change was to simply register the prototypes to the reference images by using the fat instead of the water images. The parameters used were the same as for the Demon algorithm with the water images as input which can be seen in table 3.2.

The second change was to combine the water and fat images in the registrations. This was done by first using the fat images on scale 3, 2 and 1 followed by the water images on scale 1 and 0. A list of the number of iterations done on each scale are given by table 3.3. The remaining parameters remained the same as in table 3.2.

(45)

3.3 Region definitions 35 Table 3.3:Parameters used for the Demon algorithm when both the fat and water images were used in the registrations.

Parameter Value

Start scale fat 3

End scale fat 1

Iterations per scale fat [50 50 50] Start scale water 1 End scale water 0 Iterations per scale water [50 25]

3.3 Region definitions

In order to understand the method by which the prototypes are registered to a reference image and how this leads to a segmentation of the reference image, it is important to describe which regions that are of interest to segment and how they are defined. Why these regions are of particular interest and why the mus-cles are grouped in a particular way are not relevant for this thesis, however the definitions of these regions are. The definitions are given below and images that illustrates these regions can be seen in figure 3.2. Note that there is a difference between VAT/ASAT measurements and VAT/ASAT regions. The measurements determines the total fat volume within a particular region.

VAT region

The region that determines where the visceral adipose tissue is located is defined as the region inside the peritoneum where the upper boundary is the diaphragm and the lower boundary is determined by the hip bone and sacrum.

ASAT region

The region that determines where the abdominal subcutaneous adipose tissue is located is defined by all the tissue outside the abdominal muscles where the up-per boundary is the 9th thoracic vertebrae and the lower boundary is the top of the femural head.

Left/Right Upper Leg Front

These regions, LULF and RULF, are defined as the union between the muscles sartorius, tensor fasciae latae and quadriceps femoris in the left and right leg.

(46)

Left/Right Upper Leg Back

These regions, LULB and RULB, are defined as the union between the muscles gluteus minimus, medius and maximus, semimenbranosus, semitendinosus and biceps femoris in the left and right leg.

Figure 3.2:The image to the left illustrate the upper leg masks. Blue pixels represent left upper leg back, yellow represents left upper leg front, magenta represents right upper leg back and cyan represents right upper leg front. The image to the right illustrate the VAT and ASAT regions where red pixels represent VAT and blue pixels represent ASAT. Image source: AMRA.

3.4 Dataset

This section describes the datasets used as template and reference images in this thesis. From now on, the template images is referred to as prototypes since that is the name convention in the context of full body image registration.

The prototypes and reference images were taken from the UK bio-bank which is a major international health resource. These images were obtained from an MRI scan and as mentioned in section 1.1, an MRI scan generates both a water and a fat image. The prototypes and reference images could either refer to the water or the fat images depending on the context.

A set of 29 prototypes was used to register each reference image. This was needed in order to represent many different body types and to give a more accurate seg-mentation. The prototype set was fixed to ensure a fair comparison between the algorithms and to have one parameter less to tune. The set of reference images used for evaluation consisted of 67 images with varying body composition.

(47)

3.4 Dataset 37

The size of the images are the same along the x and y axis while the size along the z axis differs. The varying sizes and the resolution of the images can be seen in table 3.4.

Table 3.4:Size and resolution of UK bio-bank images Property UK bio-bank images

x-size 224 voxels y-size 174 voxels z-size 320 - 370 voxels x-resolution 2.232 mm y-resolution 2.232 mm z-resolution 2.996 mm

Even though datasets with images of larger sizes and higher resolution were avail-able, the lack of computational power made it very time consuming to evaluate the algorithms on those datasets.

3.4.1 Region masks

Each prototype and reference image has a set of six binary masks pertaining to them that indicate which region each voxel belongs to. The binary masks corre-spond to the regions described in section 3.3. An example of a prototype and a pertaining binary mask can be seen in figure 3.3.

Figure 3.3:The image to the left illustrate the water image from a prototype and the image to the right illustrate the binary mask that indicate the left upper leg front region. Image source: AMRA.

(48)

The region masks pertaining to the prototypes are accurate and well segmented. It is however paramount to mention the precision of the ground truth masks pertaining to the reference images and their flaws. The procedure of obtaining a ground truth mask is divided into two steps:

1. The Morphon algorithm autonomously generates suggestions, i.e. binary masks for each region.

2. The suggestions are manually inspected and corrected (if needed).

Since AMRA sell body composition measurements, the suggested masks are only corrected if they contain regions that will add error to the measurement. As an example, if the ASAT mask leaks into the VAT region (which contain a fat signal), the mask will be corrected. However, if the ASAT mask leaks outside of the body (which contains no fat signal) the mask is not corrected. An image which illus-trate this effect for the muscle regions can be seen in figure 3.4. Note that the ground truth masks are corrected where the image contains a signal that will add to the muscle volume.

With this in mind, an evaluation metric which only considers the overlap between the estimated region mask generated by an algorithm and the ground truth mask pertaining to the reference image will not only be inaccurate, but also biased toward the Morphon. How this was avoided is explained in section 3.6.

Figure 3.4: An example of a flawed ground truth masks. The images to the left illustrate the front and back leg ground truth masks. The images to the right illustrate the masks after deformation with the Morphon algorithm. Blue pixels represent left upper leg back, yellow represents left upper leg front, magenta represents right upper leg back and cyan represents right upper leg front. Image source: AMRA.

(49)

3.5 Deformation and segmentation 39

3.5 Deformation and segmentation

This section explains the method by which the prototypes are registered to a reference image and how this leads to a segmentation of the reference image. To make the content easier to digest, a flowchart which describes the different steps of the method can be seen in figure 3.5.

Figure 3.5: A flowchart that describes the different steps in segmenting a reference image.

(50)

3.5.1 Prototype deformations

In order to segment a reference image, the first step was to register each proto-type to the reference image. This means that, for each protoproto-type, a displacement field was found in order to map the prototype into the reference image. The field was obtained from one of the algorithms explained in section 3.1. An example of a prototype before and after it is deformed into a reference image can be seen in figure 3.6.

Figure 3.6:This figure illustrates the effects of an image registration with the fat images in the coronal plane of a body. The image to the left is a prototype, the image in the middle is a deformed prototype and the image to the right is a reference image. Note that the prototype is more similar to the reference image after the deformation than before. Image source: AMRA.

3.5.2 Prototype selection

After the mapping between each prototype and a reference image is obtained, only a subset of prototypes is allowed to impact the segmentation. The idea is to determine which prototypes, after deformation, that are most similar to the reference image in each of the regions described in section 3.3. This is done by: First, finding the upper and lower boundary of the different regions in the refer-ence image. Second, computing the normalized cross correlation (2.18), for each region, between the deformed prototypes and the reference image. Third, the prototypes with the highest NCC values are selected to create a probability field for each region. Note that the selected prototypes might differ from region to region.

In order to find the regions of interest in the reference image, the displacement fields obtained by the registrations are applied to the region masks pertaining to each prototype. Then, the min and max z values of the deformed masks are used to determine the upper and lower boundary of these regions. The normalized cross correlation are then computed between the deformed prototypes and the reference image in the area determined by these boundaries.

(51)

3.5 Deformation and segmentation 41

Muscle masks pertaining to a prototype before and after a registration on top of reference image can be seen in figure 3.7.

Figure 3.7:These images illustrate the front and back leg masks pertaining to a prototype before and after deformation on top of a reference image. The images to the left illustrate the masks before the deformation and the images to the right illustrate the masks after the deformation. Note how the masks fit the reference image much better after the deformation. Blue pixels represent left upper leg back, yellow represents left upper leg front, magenta represents right upper leg back and cyan represents right upper leg front. Image source: AMRA.

3.5.3 Probability field

In order to determine whether a voxel in the reference image is part of a partic-ular region, a probability field P Fregionis created from the deformed masks

per-taining to the selected prototypes with the highest NCC values. The probability that a specific voxel belongs to a certain region is given by

P Fregion(v) = 1 Nregion Nregion X p=1 Mregion,p(v) (3.7)

where Mregion,pis a deformed mask for a specific region of the pthprototype, v

is a particular voxel and Nregionis the number of prototypes used to create the

(52)

3.5.4 Segmentation

The final segmentations, i.e. creating the binary masks of each region Mregionin

the reference image, are obtained by thresholding the probability fields,

Mregion(v) = (P Fregion(v) > Tregion), (3.8)

where Tregionis the threshold for a specific region.

In this thesis, the optimal combination of Nregionand Tregionwas found with an

semi-exhaustive search for the parameter that resulted in the highest mean pre-cision and recall. The term semi-exhaustive search is used since not all values of Nregionwere evaluated, only {7, 9, ..., 23}. All thresholds however, for each value

of Nregion, were evaluated.

A list of the optimal numbers of prototypes used to create the probability fields and the optimal thresholds used to segment the probability fields, for each algo-rithm and region, can be seen in table 4.1 and 4.5.

3.5.5 Disjoint regions

After the regions are segmented, there is a possibility that the binary masks over-lap. This overlap has to be dealt with in order to create disjoint regions and compute accurate body composition measurements. If a voxel is a part of two re-gions, the probability field with the highest probability, for that particular voxel, determines which region it will ultimately belong to.

3.6 Evaluation

After the regions of interest had been segmented, the registration was evaluated and compared to the ground truth. Two types of evaluation measurements, vol-ume difference and the Bland-Altman plot, were included to evaluate how the registration affects the body composition measurements. The other evaluation measurements, weighted precision and recall, were used to evaluate how the reg-istration affects the segmentation.

3.6.1 Body composition measurements

The body composition measurements used were VAT, ASAT and fat free muscle volume of the upper leg regions, all of which are described in section 2.5.1 but repeated in this section with a slightly different notation. The estimated body composition measurements are given by

(53)

3.6 Evaluation 43

V ATest=

X

∀_voxels

f · MV AT ,est· Mbody· Vvoxel, (3.9)

ASATest=

X

∀_voxels

f · MASAT ,est· Mbody· Vvoxel, (3.10)

LU LFest=

X

∀_voxels

(f < 0.5) · MLU LF,est· Mbody· Vvoxel, (3.11)

RU LFest=

X

∀_voxels

(f < 0.5) · MRU LF,est· Mbody· Vvoxel, (3.12)

LU LBest=

X

∀_voxels

(f < 0.5) · MLU LB,est· Mbody· Vvoxel, (3.13)

RU LBest=

X

∀_voxels

(f < 0.5) · MRU LB,est· Mbody· Vvoxel. (3.14)

The fat-concentration, f , is the intensity in a particular voxel of the reference im-age. MV AT ,est, MASAT ,est, MLU LF,est, MRU LF,est, MLU LB,est and MRU LB,estare the

estimated binary masks generated by the method explained in section 3.5 and correspond to (3.8). Mbody is a binary body mask of the reference image. Note

that each voxel within any of the muscle regions with a fat concentration less than 0.5 is considered to be a fat-free muscle voxel.

To determine how accurately each algorithm computed these measurements, vol-ume differences and Bland-Altman plots were used. These measurements tell how accurately the algorithm estimate the correct volume, but not necessarily how well the algorithm determine the correct segmentation.

The total volume difference vol,tot is given by the difference between the

esti-mation of the volume and the ground truth volume in a particular region. The equation is given by

vol,tot = bcmest−bcmgt, (3.15)

where bcmestis the estimated measurement and bcmgt is the ground truth

mea-surement. Bcmestcan refer to either V ATest, ASATest, LU LFest, RU LFest, LU LBest

or RU LBestgiven by (3.9 - 3.14). Bcmgtis generated by the same equations where

the estimated binary mask has been substituted by the ground truth binary mask for each region.

The Bland-Altman plot explained in section 2.5.1 was used to give a visual rep-resentation of the difference between an estimated and a ground truth body com-position measurement for each reference image in the evaluation set.

Evaluation of 3D MRI Image Registration Methods

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2017