Nonrigid Registration of Dynamic Contrast-Enhanced MRI ... - TSpace

Nonrigid Registration of Dynamic Contrast-Enhanced MRI data using Motion Informed Intensity Corrections

by

Anthony Lausch

A thesis submitted in conformity with the requirements for the degree of Master of Science Graduate Department of Medical Biophysics University of Toronto

Copyright

© 2011 by Anthony Lausch

Abstract Nonrigid Registration of Dynamic Contrast-Enhanced MRI data using Motion Informed Intensity Corrections Anthony Lausch Master of Science Graduate Department of Medical Biophysics University of Toronto 2011 Effective early detection and monitoring of patient response to cancer therapy is important for improved patient outcomes, avoiding unnecessary procedures and their associated toxicities, as well as the development of new therapies. Dynamic contrast-enhanced magnetic resonance imaging shows promise as a way to evaluate tumour vasculature and assess the efficacy of new anti-angiogenic drugs. However, unavoidable patient motion can decrease the accuracy of subsequent analyses rendering the data unusable. Motion correction algorithms are challenging to develop for contrast-enhanced data since intensity changes due to contrast-enhancement and patient motion must somehow be differentiated from one another. A novel method is presented that employs a motion-informed intensity correction in order to facilitate the registration of contrast-enhanced data. The intensity correction simulates the presence or absence of contrast agent in the image volumes to be registered in an attempt to emulate the level of contrast-enhancement present in a single reference image volume.

ii

Acknowledgements I would like to take this opportunity to thank everyone who made this project possible. First, I owe a great deal of gratitude to my supervisor Dr. Anne Martel for taking me on as a graduate student and providing me with the opportunity to work on such an interesting project. Her adept supervision combined with her personable and laid-back style created a perfect environment in which to learn and grow as a researcher. I could not have asked for a better supervisor. Similarly, I would like to thank my committee members, Drs. Jason Lerch and Kristy Brock for their very helpful perspectives and guidance throughout the course of the project. The members of the Martel research group - Cristina, Hatef, Danoush, Rushin, Jacob, Gary, James and Mehran - provided me with invaluable support and constructive feedback. Special thanks to Mehran for all of our fun discussions on registration and the underlying mathematics. Colleen Bailey was extraordinarily helpful in providing me with and explaining the tools I needed to test my algorithm from a clinical perspective. Her continually precise efforts on my behalf significantly contributed to the timely completion of my project. Dr. Greg Stanisz and Dr. Georg Bjarnason also contributed in this fashion by providing me with the much needed datasets to test my algorithms. Rafal Janik was also very helpful and always available to answer questions regarding MRI. Finally, a big thank-you to all of my friends and family for general support throughout the project. I am extremely fortunate to have such a great group of people cheering me on and keeping me balanced.

iii

Contents

1 Introduction

1

1.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2

DCE-MRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2.2

Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.3

. . . . . . . . . . . . . . . . . . . . . . . . .

8

1.3.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.3.2

Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.3.3

Other Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

1.4

Motion Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

1.5

Image Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

1.5.1

Distance Measures . . . . . . . . . . . . . . . . . . . . . . . . . .

16

1.5.2

Regularizers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

1.5.3

Transformation Models . . . . . . . . . . . . . . . . . . . . . . . .

20

1.5.4

Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

Registration of Contrast-Enhanced Data . . . . . . . . . . . . . . . . . .

25

1.6.1

Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

1.6.2

Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

Hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

1.6

1.7

Pharmacokinetic Modelling

iv

2 Registration Using Intensity Corrections

30

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

2.2

Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

2.2.1

Patient Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

2.2.2

Image Registration Framework

. . . . . . . . . . . . . . . . . . .

32

2.2.3

Reference Selection . . . . . . . . . . . . . . . . . . . . . . . . . .

33

2.2.4

Intensity Correction

. . . . . . . . . . . . . . . . . . . . . . . . .

33

2.2.5

Estimating Mmax . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

2.3.1

Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

2.3.2

Pharmacokinetic Modelling . . . . . . . . . . . . . . . . . . . . .

41

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

2.4.1

Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

2.4.2

Pharmacokinetic Modelling . . . . . . . . . . . . . . . . . . . . .

45

2.5

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

2.6

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

2.3

2.4

3 Conclusions and Future Work

55

3.1

Significance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

3.2

Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

3.2.1

Blurring Kernel Definition and Motion Estimation . . . . . . . . .

57

3.2.2

Resolution and Registration Accuracy . . . . . . . . . . . . . . . .

59

3.2.3

Application and Testing . . . . . . . . . . . . . . . . . . . . . . .

60

3.2.4

Effects of Registration and Motion on Analysis . . . . . . . . . . .

60

Summary of Thesis Contributions . . . . . . . . . . . . . . . . . . . . . .

61

3.3

Appendix

62

A Pharmacokinetic Model Fitting

64 v

B Principal Components Analysis and Imaging

67

Bibliography

69

vi

List of Figures 1.1

Breast DCE-MRI subtraction image with a potential malignancy. . . . .

5

1.2

Typical DCE-MRI intensity-time curve. . . . . . . . . . . . . . . . . . . .

6

1.3

Typical 2-dimensional SPGR pulse sequence.

. . . . . . . . . . . . . . .

7

1.4

Two-compartment Tofts model. . . . . . . . . . . . . . . . . . . . . . . .

9

1.5

Tofts model fitting procedure for DCE-MRI data. . . . . . . . . . . . . .

11

1.6

Effect of motion on ROI tissue content. . . . . . . . . . . . . . . . . . . .

15

1.7

Effects of motion on mean value ROI intensity-time curves. . . . . . . . .

15

1.8

Multi-resolution representation of an abdominal DCE-MRI image. . . . .

24

1.9

Mock image of a kidney-like organ with partial enhancement after time ∆t. 26

2.1

Intensity correction for a simulated enhancing necrotic tumour. . . . . . .

35

2.2

Simulated enhancing tumour images and intensity correction. . . . . . . .

36

2.3

Quantitative validation pipeline . . . . . . . . . . . . . . . . . . . . . . .

40

2.4

Distance of voxels from a reference image volume before and after registration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

2.5

Relative changes in chi-squared scores for model fits due to registration. .

45

2.6

Chi-squared scores after registration for voxel-based fits. . . . . . . . . .

46

2.7 K trans parameter maps before and after registration. . . . . . . . . . . .

47

2.8 ve parameter maps before and after registration. . . . . . . . . . . . . . .

48

2.9

49

Intensity-time curve before and after registration. . . . . . . . . . . . . . vii

List of Abbreviations AIF Arterial Input Function AP Anterior-Posterior CT Computed Tomography DCE Dynamic Contrast Enhanced EES

Extracellular Extravascular Space

FWHM Full Width at Half-Maximum jPDF

Joint Probability Distribution Function

l-BFGS Limited memory Broyden-Fletcher-Goldfarb-Shanno LR Left-Right MI Mutual Information MRI

Magnetic Resonance Imaging

NCC Normalized Cross Correlation NGF Normalized Gradient Fields NMI

Normalized Mutual Information

PCA

Principal Components Analysis

PET

Positron Emission Tomography

R RF

Reference Image Radiofrequency

ROI Region of Interest SI Superior-Inferior SPGR Spoiled Gradient Echo SSD Sum of Squared Differences T

Moved Image

TE Echo Time TR Repetition Time US Ultrasound viii

Chapter 1 Introduction

1

Chapter 1. Introduction

1.1

2

Overview

Cancer is a leading cause of death world-wide, accounting for roughly 13% of all deaths [63]. In Canada, cancer accounts for nearly 30% of all deaths and is the most significant contributor to the number of premature years of life lost [14]. Common cancer treatment approaches generally involve chemotherapy, radiotherapy, surgery, or potentially a combination of these treatments. In-vivo imaging has become an essential tool for the noninvasive assessment of patient response to treatment. Anatomical and functional imaging can be performed to assess the current size and number of tumours within the body. Computed tomography (CT) and magnetic resonance imaging (MRI) are particularly useful for identification of tumours and measurements of tumour size, while functional imaging modalities such as nuclear medicine and positron emission tomography (PET) scans can reveal metastases throughout the body. Early detection and monitoring of response to therapy is important for many reasons. Foremost amongst these is the potential to modify therapy in the event of partial response or disease progression. Ideally this should occur as soon as possible in order to avoid unnecessary treatments and their associated toxicities as well as to improve outcomes under modified treatment regiments. Additionally, effective methods for early evaluation of response are helpful for the development of new treatments. Large and expensive clinical trials are generally required to verify efficacy of new therapies. Enhanced evaluation capabilities may reduce the number of participants required to gain useful statistics on patient response thus decreasing the time and expense required to develop new therapies. Consequently there has been a large push to develop new techniques to evaluate and monitor response to treatment [36]. Changes in gross measurements of tumour size have been correlated with patient outcome [20], however these changes can take time to manifest or may not occur at all depending on the type of treatment. Characterizing changes in tumour vasculature caused by treatment may offer a promising new alternative to evaluate and monitor


3

patient response earlier on in treatment. Local vasculature plays an essential role in tumour growth, metabolism, and tissue repair [31, 43]. As a result, vasculature is often affected by or essential to treatment. For example, some chemotherapy drugs seek to inhibit tumour growth by disrupting the formation of new tumour vasculature (antiangiogenics)[43]. Alternatively, radiotherapy may be less effective when applied to poorly perfused tumours with limited oxygen supply [61]. Dynamic contrast-enhanced (DCE) imaging modalities such as DCE-MRI, DCE-CT, and DCE ultrasound (DCE-US) show promise as potential methods to evaluate the state of tumour vasculature [48, 30, 32, 55]. With these approaches, the patient is repeatedly imaged as an injected contrast-agent flows through the patient vasculature. The contrastagent creates signal intensity changes that are visible in the acquired images. Depending on the modality, the contrast-agent dynamics observed in the images can then be linked to properties of the vasculature and surrounding tissue such as vessel permeability and blood flow [36, 30, 55]. Changes in these quantities may be helpful for characterizing response. Since imaging occurs repeatedly over time and the patient will move (e.g. breath), tissue in images acquired at one point in time will not necessarily be located in the same position in images acquired at another point in time. In order to accurately characterize a tissue of interest, the tissue’s position must be known at all times. Image registration is a post-processing technique that can be used to remove motion from dynamic image datasets. The goal of image registration is to find a transformation that will move the structures in one image such that they align spatially with the same structures contained in a different image. In effect, image registration tracks the position of structures in the images and then outputs a modified set of images where the structures have been moved back to some original reference location. Although image registration has been developed and successfully applied within many imaging contexts, it is challenging to develop for contrast-enhanced image data. Struc-


4

tures in contrast-enhanced datasets both move and change in appearance which significantly increases the difficulty in automatically establishing structural correspondences in different images. In this thesis, a novel method for registration of contrast-enhanced images is presented and tested using abdominal DCE-MRI data. Quantitative accuracy is demonstrated using a realistic simulation and potential clinical utility is investigated through looking at the effects of registration on DCE-MRI analyses. In the proceeding sections of this chapter, DCE-MRI and related analysis methods are explored in further detail. The problem of patient motion is then illustrated within this context. Common image registration techniques are described in depth, followed by a discussion of why these methods are difficult to implement for contrast-enhanced data. In Chapter 2, the novel registration method is introduced, tested, and discussed. Chapter 3 concludes the document and suggests future research directions.

1.2 1.2.1

DCE-MRI Overview

Dynamic-contrast enhanced magnetic resonance imaging (DCE-MRI) involves the acquisition of multiple MR images as an injected contrast-agent flows through the body. Paramagnetic contrast agents produce signal increases in T1-weighted images by shortening the T1-relaxation times of neighbouring hydrogen atoms. The dynamics of the contrast agent and the resulting local signal enhancement are related to several different factors that characterize local vascular structure and function. The concentration of the contrast agent in the arterial blood supply, tissue perfusion, capillary surface area and permeability, and the volume of extracellular extravascular space (EES) all influence the contrast agent dynamics observed with DCE-MRI [48]. Consequently, DCE-MRI has been used to investigate the state of tumour vasculature. This technique has been used


5

in many different sites including, but not limited to, breast, kidney, liver, prostate, brain, lung, heart and bone [45, 66, 47, 68, 42, 53, 24, 2, 50, 49, 34, 15, 6]. Unlike DCE-CT, DCE-MRI is non-ionizing and can be performed repeatedly without worrying about accumulated radiation dose. It also provides excellent soft tissue contrast compared to DCE-CT and DCE-US, but can suffer from lower temporal resolution. In order to grow, tumours must have access to increased blood supply to meet their increasing metabolic needs (e.g. oxygen, nutrients, waste removal). Cancer cells release factors (e.g. vascular endothelial growth factor) that stimulate the formation of new blood vessels which can provide for these requirements and support further growth [31]. In tumours, the formation of new blood vessels (angiogenesis) is highly disorganized and results in abnormal leaky vasculature. These differences in vasculature influence local enhancement dynamics observed in DCE-MRI. An image taken prior to the injection of the contrast agent can be subtracted from the post-injection images in order to isolate the contrast-agent dynamics. Suspicious regions of interest (ROIs) can be identified based on visual enhancement characteristics and further investigated for potential malignancy. Breast DCE-MRI provides a strong illustrative example (Figure 1.1).

Figure 1.1: Breast DCE-MRI images. The post-injection image appears on the left and the subtraction image is on the right. A potential malignancy is circled in white. A more quantitative analysis can be performed by investigating the signal intensity-time curves of suspicious ROIs (Figure 1.2). In some sites, such as the breast, the shape of the


6

intensity-time curve can give insight into whether a tumour may be benign or malignant [35].

Figure 1.2: Typical intensity-time curve. I) Pre-injection phase II) Contrast agent washin phase III) Contrast agent wash-out phase. Using the intensity-time curve and knowledge of the MR pulse sequence, the concentration of the contrast agent can be calculated as a function of time. A pharmacokinetic model can then be fit to the concentration-time curve. The fitting procedure yields parameters that relate to vascular permeability, blood volume, and the EES. New anticancer drugs which seek to inhibit tumour growth by disrupting angiogenesis affect local vasculature. Consequently DCE-MRI shows promise as a new method to evaluate the efficacy of these drugs [36, 48]. Standard techniques to evaluate therapeutic response generally rely upon investigating changes in gross measurements of tumour size [20]. However, tumour size changes due to reduced blood supply can take time to manifest in anatomical images. In addition, antiangiogenic drugs may only stop further tumour growth rather than actually reduce tumour size. In contrast, DCE-MRI provides a direct method for evaluating the mechanism of action of these drugs and could potentially be used to monitor patient response. Effective evaluation can lead to treatment course changes and improved patient outcomes.

7


1.2.2

Techniques

T1-weighted DCE-MRI data can be acquired using a variety of different pulse sequences. Many of these correspond to modifications of a spoiled gradient echo sequence (SPGR, Figure 1.3). First, an α° radiofrequency (RF) pulse introduces a transverse component to the magnetization of the spins. The slice-select gradient (Gslice ) is simultaneously turned on, along with a varying phase-encode gradient (Gphase ) to permit sampling of 2dimensional k-space. The readout gradient (Greadout ) is also turned on in order to dephase the spins. At a time TE (echo time) later, a Greadout pulse with opposite polarity is applied to bring the spins back into phase. At this time the signal is recorded. A pseudorandom phase (αphase ) is added to the RF pulses upon each repetition to produce spatially independent phase shifts.

Figure 1.3: Typical 2-dimensional SPGR pulse sequence. Smaller flip angles can permit shorter repetition times (TR) that in turn facilitate faster image acquisition. Fast image acquisition can be important for characterizing contrast-agent dynamics. For short TR, a spoiler gradient is required in order to avoid the contamination of subsequent repetitions with residual transverse magnetization. The SPGR signal equation is given by, M0 (1 − e−( T 1 ) )sin(α) TR

s(t) =

1 − e−( T 1 ) cos(α) TR

(1.1)


8

where s is the recorded signal, TR is the pulse sequence repetition time, T1 is the longitudinal relaxation time (spin-lattice), M0 is the equilibrium longitudinal magnetization, and α is the flip angle. A similar pulse sequence can be used for 3-dimensional SPGR imaging by using a Greadout gradient with a different magntiude upon each repetition. Though small flip angles and shorter TR can permit faster imaging, an inherent trade-off exists between temporal and spatial resolution. As α and TR are reduced the resulting signal decreases which contributes to lower spatial resolution and contrast. To improve signal and contrast within the context of extremely small flip angles(< 5 °) and TR (around 1s), an initial 180° RF preparation pulse can be applied prior to the SPGR pulse sequence. This technique is particularly useful within contexts where extremely rapid acquisition is required (e.g. cardiac). A drawback to this approach is the limited field of view that can be scanned due to the extremely short TR time.

1.3 1.3.1

Pharmacokinetic Modelling Overview

Within the context of contrast-enhanced imaging, pharamcokinetic modelling endeavours to characterize the dynamics of the contrast agent as it flows through the body. By fitting a model to the concentration-time curves of ROIs within the image data, information pertaining to blood flow and vascular permeability can be acquired [48]. Although pharmacokinetic modelling is not the focus of this thesis, a general overview is needed in order to understand how registration may improve analysis of clinical data. Proceeding discussions focus on a single model for illustrative purposes, since other models are all affected by motion in a similar fashion. Tofts et al [59] introduced what has become widely known as the standard twocompartment Tofts model. They summarized many of the different pharmacokinetic

9


models that had previously been developed and derived a common set of characteristic parameters and naming conventions. The model assumes that there are two compartments that correspond to the vasculature (blood plasma) and the EES (tissue)(Figure 1.4).

Figure 1.4: Two-compartment Tofts model. The contrast agent diffuses from the blood plasma to the tissue and from the tissue back to the blood plasma at rates described by the volume transfer constants K trans and kep respectively. The arterial input function (AIF) specifies the concentration of the contrast agent in the blood plasma at some initial time point after the injection of the contrast agent. The concentration of the contrast agent in blood plasma (Cp ) and in the tissue (Ct ) is then modeled by an ordinary differential equation given by dCt = K trans Cp − kep Ct dt

(1.2)

dCt Ct = K trans (Cp − ) dt ve

(1.3)

or equivalently,

where ve correspond to the fraction of the tissue volume occupied by the EES and kep = K trans /ve . The AIF serves as an initial condition which is used to solve the differential

10


equation. Generally it is assumed the AIF is a short arterial pulse (delta function) in contrast-agent concentration that brings the concentration in the vasculature from zero to some non zero value Cp0 , at time t = 0. The solution to equation 1.2 for a delta function AIF is given by Ct =

K trans Cp0 (1 − e−kep t ). kep

(1.4)

A small modification of the standard Tofts model can improve the modelling of the concentration of the contrast agent in the plasma compartment. Instead of assuming that the concentration is described by an impulse function that increases Cp at t = 0 from 0 to Cp0 , the actual time-varying concentration can be measured directly in a nearby artery at all time points. Cp (t) is then explicitly defined over the duration of the scan, rather than by an impulse function that is scaled by some initial peak estimate of Cp (t = 0). The interpretation of the parameter K trans can vary depending on the conditions of vascular permeability that are assumed. For high permeability K trans is related to the perfusion of whole blood per unit mass of tissue [33], while for low permeability it is related to the permeability surface area product of the vasculature per unit mass of tissue [59]. For a mix between high and low permeability conditions, K trans is related to both perfusion and the fractional reduction of the contrast agent concentration as it passes through the capillaries [59].

1.3.2

Fitting

The fitting of a pharmacokinetic model to DCE-MRI data involves multiple steps which are summarized in figure 1.5 for a Tofts model with a measured AIF. In step 1, the concentration of the contrast agent in the blood plasma, Cp (t), is determined in order to define the AIF for all time. First, the time-varying average pixel (voxel) intensity is calculated for an ROI placed within an artery that feeds the tissue of interest. These values serve as the SPGR signal for the blood plasma, Sp (t). Knowing Sp (t), the flip angle α, the repetition time TR, and the equilibrium longitudinal

11


Figure 1.5: Tofts model fitting procedure for DCE-MRI data. The concentration of the contrast agent is denoted by [CA]. relaxation magnetization Mp0 , the SPGR equation can be solved for the T1 of the blood plasma, T 1p (t). Mp0 must be previously determined by performing some preliminary MRI experiments (further discussed in Appendix A). After T 1p (t) is determined, Cp (t) is then given by Cp (t) =

1/T 1p (t) − 1/T 1p0 r(1 − Hct)

(1.5)

where T 1p0 is the relaxation time of the plasma when the contrast agent is not present. r is the relaxivity of the contrast-agent which characterizes its ability to change the T1 relaxation times of neighbouring spins. If the ROI used to determine Cp (t) is not directly within the tissue of interest, a small time offset can be applied to Cp (t) to account for the time it takes for the contrast agent to flow from the ROI to the tissue. In step 2, initial values of the pharmacokinetic model parameters are guessed or esti-

12


mated. For the standard Tofts model these parameters would correspond to K trans and ve . These values serve as a starting point for subsequent model fitting and optimization procedures. In step 3, the estimated parameters are used along with the pharmacokinetic model to predict the concentration of the contrast agent in the tissue of interest. In the case of the standard Tofts model with a delta function AIF, equation 1.2 can be solved explicitly as in equation 1.4. However, for a measured time-varying AIF, the model must be integrated numerically to compute Ct (t). In step 4, the SPGR signal that would result from the predicted concentration of contrast agent in the tissue is calculated (Stpred ). In order to do this, T 1t (t) is calculated from the predicted Ct (t) and then substituted into the SPGR equation. In step 5, the actual SPGR signal from the tissue of interest in the DCE-MRI data (St ) is compared to Stpred (t). A common approach is to use the chi-squared or reduced chi-squared score as a metric: χ2 =

∑ (St − Stpred )2 i

χ2red =

σ2 1 χ2 DOF

(1.6)

(1.7)

where DOF is the degrees of freedom and is equal to the number of observations less the number of fitting parameters less 1. St is sampled in different ways depending on whether mean-value ROI or voxel-based analysis is being performed. Consider the task of modelling the vascular characteristics of a spherical tumour using the standard Tofts model [59]. For mean-value ROI analysis, St corresponds to the average intensity of the voxels in the tumour ROI at each time point. After the fitting procedure is finished, this approach yields a single average estimate of K trans and ve for the entire tumour. For voxel-based analysis, St corresponds to the time-varying intensity of a single voxel within the tumour ROI. After fitting, this yields a single estimate of K trans and ve for the voxel. The voxel-based fitting procedure is then repeated for all of the voxels in the tumour


13

ROI resulting in unique estimates of K trans and ve for each voxel. Finally in step 6, the value of the metric is optimized by updating the estimates of the model parameters in an iterative optimization loop. Minimization of the chi-squared scores increases the agreement between the observed and model predicted SPGR signal for the tissue. If the model accurately predicts the true SPGR signal, then the associated model parameters are also assumed to be accurate (e.g. K trans , ve ).

1.3.3

Other Approaches

The standardized Tofts model and its variations have enjoyed broad use due to their simplicity and fitting parameters that have been shown to have some diagnostic and prognostic utility [65]. However, many different pharmacokinetic models have been developed and tested on DCE-MRI data [65]. Brix et al. [9] introduced the two-compartment exchange model which is an extension of the Tofts model. The relevant fitting parameters include the blood plasma compartment fractional volume (PV), the EES compartment fractional volume (EV = ve ), the plasma flow in the vasculature (PF), and the EES flow (EF). It is assumed that the rate of flow into and out of the EES is equal to EF at all times. The addition of PF and PV as fitting parameters permits more realistic contrast agent dynamics. With the Tofts model, the concentration of the contrast agent in the plasma is equivalent to the AIF at all times [38] whereas the Brix variation models this concentration more explicitly. For example, the Brix model can permit widening of the contrast-agent bolus as it passes through the vasculature [38]. Simpler, semi-quantitative models have also been investigated. Some examples include measurements of the area under the intensity-time curve, time-to-peak enhancement, and slope of the washout portion of the intensity-time curve [65]. Although these measurements do not directly reflect vascular parameters, they are easy to implement and are correlated with a mix of parameters such as blood volume and permeability [65].


1.4

14

Motion Problem

In practice there is always some amount of motion present in DCE-MRI data since the dynamics of the contrast-agent must typically be characterized over a period of minutes. Various sources of patient motion exist on such time scales. The heart contracts and deforms as it pumps blood and the diaphragm moves and displaces surrounding tissue and organs during respiration. Respiration induced motion of abdominal organs such as the liver, kidneys and spleen has been observed to be anywhere between 1 mm and 50 mm depending on the individual patient [64]. Patients may also become uncomfortable in the scanner and shift to reposition themselves. Conversely, they may become more comfortable as they become accustomed to the procedure leading to involuntary muscle relaxation (e.g. pectoral muscles). Consequently motion is present within many of the contexts in which DCE-MRI is utilized. Motion can significantly reduce the accuracy of DCE-MRI analysis and can often render image data unusable [48]. DCE-MRI data have frequently been omitted from clinical trials due to motion artifacts [48]. The primary obstacle relates to the definition of the ROIs used for pharmacokinetic modelling. An ROI delineated on a set of DCEMRI images at one point in time does not necessarily contain all of the same tissue at another point in time due to patient motion (Figure 1.6). Consequently, the resulting intensity-time curves may not reflect the precise contrast-agent dynamics of the tissue under investigation (Figure 1.7). This problem can be partially addressed through mean-value ROI analysis, which can be robust in the presence of motion but suffers from insensitivity when the underlying tissue is heterogeneously enhancing [36]. In addition, this approach becomes less robust as the size of the ROI approaches the magnitude of the motion (Figure 1.7). If the two are equivalent, then at some points in time, the underlying tissue of interest will be entirely outside of the defined ROI and replaced by some other tissue with potentially different enhancement patterns. A voxel-by-voxel analysis can more accurately reflect the


15

Figure 1.6: Coronal DCE-MRI images. Renal tumour delineated in red (left). Displaced renal tumour with previous tumour delineation superimposed (right).

Figure 1.7: Effects of motion on mean value ROI intensity-time curves. The intensitytime curve for a renal tumour undergoing small magnitude motion relative to the tumour size is shown on the left. The intensity-time curve for a renal tumour undergoing larger magnitude motion relative to the tumour size is shown on the right. heterogeneity of ROIs, but is even more sensitive to motion as the voxels are in effect, very small ROIs. Abdominal DCE-MRI is particularly affected by respiratory motion in the superiorinferior (SI) and anterior-posterior (AP) directions. SI motion of the kidneys, liver/pancreas, and spleen can exceed 15 mm, 20 mm, and 13 mm respectively [13, 8]. Less motion is typically observed in the AP directions (5-9 mm) [13, 8]. Abdominal motion is also nonrigid resulting in deformation of structures.

16


1.5

Image Registration

Image registration is an image processing technique that can be used to correct for motion within or between image datasets. It involves rigidly or non-rigidly transforming images in an attempt to realign structures that have moved during the time that passes between image acquisitions. Typically a single image is selected to serve as a reference to which all other images are then registered to. If registration is accurate, then registering all of the images to the same reference automatically registers all of them to each other. There are several main components that are common to most registration algorithms. These include, 1) Distance Measures 2) Regularizers and 3) Transformation Models and 4) Optimization. For the ensuing discussion we shall consider the task of finding a transformation y that aligns some moved image T to a reference image R. The distance measure quantifies the similarity between R and T while the regularizer quantifies the suitability of the candidate transformation y. The transformation model provides a method for describing and applying the candidate transformation y. Finally, the optimization procedure searches for a transformation that optimizes the distance measure and regularizer terms with respect to the transformation model. The following concepts can all be directly extended to the registration of three dimensional image volumes.

1.5.1

Distance Measures

Distance measures quantify the similarity between the reference image R and the moved image T . A large number of intensity-based distance measures have been developed [44]. The sum of squared differences (SSD) and mutual information (MI) [16, 60, 44] distance measures are perhaps the most commonly used. The SSD is given by SSD(T, R) =

∑

(T (y(xi )) − R(xi ))2

(1.8)

i

where T (y(xi )) is the pixel or voxel intensity at position xi , of the image T that has been deformed by the candidate transformation y. R(xi ) is the pixel or voxel intensity

17


of the image R at position xi . The SSD assumes that the intensity differences between T and R are solely caused by motion and that the intensity of a given tissue type does not change. As T and R become more and more aligned, the SSD decreases with perfect alignment corresponding to a value of zero. Due to the assumption of stationary pixel intensities, the SSD is generally only applicable when registering images acquired using the same modality. MI [16, 60] does not assume that similar structures in the two images have the same pixel intensities. Consequently it is well-suited for use in multi-modal registration. The mutual information is given by M I(T, R) = H(T ) + H(R) − H(T, R) where H is the entropy H(R) = −

∑

(1.9)

pR (i)log(pR (i))

(1.10)

pT (i)log(pT (i))

(1.11)

i

H(T ) = −

∑ i

H(T, R) = −

∑∑ i

pT,R (i, j)log(pT,R (i, j)).

(1.12)

j

pR (i) and pT (i) are the probability distribution functions for the intensities in R and T evaluated at an intensity i. pR,T (i, j) is the joint probability distribution function for R and T evaluated at intentities i and j. The indices are summed over all possible intensity values in the images. It is assumed that optimal image alignment is achieved when the MI is maximized. This corresponds to when the marginal entropies of T and R are maximized and the joint entropy is minimized. However, Studholme et al [56] pointed out that in some overlapping image regions, maximization of MI in effect only corresponds to maximizing the marginal entropies. In these cases, the maximal MI, may not correspond to optimal image alignment since the joint entropy function is not minimized. Consequently they introduced the normalized mutual information (NMI) distance measure which is invariant to overlapping regions

18

Chapter 1. Introduction and changes in the marginal entropies. The NMI is given by N M I(T, R) =

H(T ) + H(R) H(T, R)

(1.13)

Maximization of the NMI implies minimization of the joint entropy function with respect to the marginal entropies and can permit registration that is invariant to image overlap. The normalized cross correlation distance measure (NCC) is another distance measure that does not assume stationary image intensities [44]. Instead the NCC assumes that there is a linear dependence between the pixel intensities in the two images to be registered. Image alignment is performed by maximizing the NCC. The NCC is given by, N CC(T, R) = where ⟨T, R⟩ =

∑ i

||T || = ||R|| =

⟨T, R⟩ ||T ||||R||

(1.14)

T (y(xi ))R(xi )

(1.15)

√ ⟨T, T ⟩

(1.16)

√

⟨R, R⟩

(1.17)

Similar to MI/NMI and NCC, the normalized gradient fields (NGF) distance measure [26] also does not assume stationary image intensities. Instead it is assumed that the edges of structures in T and R are preserved, but deform in some fashion. The NGF distance measure is given by, N GF (T, R) =

∑

1 − (n[T (y(xi ))]n[R(xi )])2

(1.18)

i

where ∇T n[T ] = √ |∇T |2 + η 2

(1.19)

The edge parameter η is a constant used to specify what constitutes an edge and what constitutes noise. Gradients with magnitudes that approach the edge parameter are suppressed as they do not convey useful physical information. Since the strength of the


19

gradients could be different in the two images to be registered, the gradient is normalized so that the distance measure is independent of gradient strength.

1.5.2

Regularizers

The regularizer quantifies the suitability of the candidate transformation y and restricts the search to physically realizable solutions. Many different transformations could yield a local minimum in the distance measure, but such transformations may not be physically realistic. For example, a transformation is clearly not suitable if it appears to break a bone or tear tissue. We also would not expect severe tissue deformation and volume changes for some sites. To address these challenges a variety of different regularizers have been developed [11, 29, 7, 22, 52]. By optimizing the regularizer at the same time as the distance measure, alignment occurs in a fashion that is restricted by the regularizer. Many regularizers focus on enforcing the smoothness of transformations in order to avoid unphysical tearing and bending phenomena. Several examples include the diffusion, curvature, and thin-plate spline regularizers [29, 22, 7]. The diffusion and curvature regularizers correspond to the sum of squares of the first order and second order partial derivatives of the transformation respectively. Consequently they penalize aggressive local changes in the ouput transformation. Thin-plate spline regularization ensures smooth transformations by quantifying the energy required to bend a thin-plate spline. This energy is described by the inner product of the Hessian of the deformation with itself. Rohlfing et al [52] proposed a localized volume preserving regularizer for registration of breast DCE-MRI data. Since tissue is generally incompressible, they assume that volume changes should be minimal. The regularizer measures the difference between the Jacobian determinant of the transformation and unity. Values over and under unity correspond to local regions of expansion and contraction respectively. Broit [11] introduced the elastic regularizer which relates the deformation of the image to the energy required to deform a homoegeneous isotropic linear elastic material. The

20

Chapter 1. Introduction elastic potential is given by S(y) =

1∑ µ < ∇y, ∇y > +(λ + µ)(∇ · y)2 2

(1.20)

where ∇y and ∇·y are the gradient and divergence of y respectively. The elastic material is characterized by the Lame parameters µ and λ. µ represents the shear modulus and λ is related to the force it takes to compress the material. Minimization of the elastic regularizer ensures that the resulting transformation treats the tissue in the images similar to a linear elastic solid described by the two Lame parameters. Although tissue can behave like a nonlinear elastic solid, this serves as a reasonable approximation, especially for smaller deformations. Consequently, the elastic regularizer is very commonly used for registration [17, 44]. It can restrict volume changes and ensure the smoothness of the transformation.

1.5.3

Transformation Models

Image registration is usually achieved by finding some transformation that minimizes a cost function (J). J is typically composed of a weighted sum of the distance measure and the regularizer. A transformation found by minimizing J, improves the alignment of the images in a way that is restricted by the regularizer. J = D(T (y), R) + αS(y)

(1.21)

Here T (y) corresponds to the image T deformed by transformation y, and α is a scalar weighting parameter. A model is needed to mathematically describe the transformations that are the focus of optimization. The model specifies how features in T can be mapped to align with features in R. Many different parametric and non-parametric models exist including, but not limited to, rigid, affine, free form b-spline, demons, optical flow, and finite element modelling [17, 54, 58, 29, 10].

21


A simple motivating example is given by the rigid deformation model. Here, T is only permitted to translate and rotate in order to align with R. In three dimensions this corresponds to a 6 parameter model where 3 parameters describe translation in the three coordinate directions,and the other 3 parameters describe rotation about the three coordinate axes. The images are then aligned by finding a set of rotation and translation parameters that minimizes the cost function. A rigid transformation model can only describe a very specific set of transformations. A less restrictive free-form model was implemented by Rueckert et al [54]. In this model, transformations are described by manipulating an underlying grid of control points. Moving a control point from one location to another results in that portion of the image moving in the same way. The control points are spanned by B-splines (polynomial basis functions) that are used to define the transformation at points that are located in between the control points. For a 3-dimensional set of control points ϕi,j,k where i ϵ [1, nx ], j ϵ [1, ny ] and k ϵ [1, nz ] the transformation f at a point (x, y, z) is given by f (x, y, z) =

3 ∑ 3 3 ∑ ∑

Bl (u)Bm (v)Bn (w)ϕi+l,j+m,k+n

(1.22)

l=0 m=0 n=0

where i = ⌊x/nx ⌋ − 1, j = ⌊y/ny ⌋ − 1, k = ⌊z/nz ⌋ − 1, u = x/nx − ⌊x/nx ⌋, v = y/ny − ⌊y/ny ⌋ and w = z/nz − ⌊z/nz ⌋. The basis functions are given by, B0 (u) = (1 − u)3 /6

(1.23)

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

(1.24)

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

(1.25)

B3 (u) = u3 /6

(1.26)

Images are aligned by finding the optimal set of control point locations (parameters) that define a transformation that minimizes the cost function. Since any number of control points can be initially defined, the number of optimization parameters is arbitrary allowing for more complicated transformations. For instance, Rueckert et al [54] considered


22

control point spacing of 10mm, 15mm, and 20mm, for registration of 3D breast DCEMRI data. However, increasing the number of control points (parameters) comes at the cost of computational efficiency. The purpose of the basis functions is to allow the use of fewer control points to describe the transformation. The basis functions can restrict the types of transformations that can be considered but become less restrictive as the number of control points increases due to their local support.

The least restrictive transformation model involves optimizing the positions of a very large grid of control points, similar to the free-form model but with far more control points and without the use of basis functions. Here, the number of control points tends to approach the number of pixels or voxels in the image to be registered. Such models are sometimes referred as ’non-parametric’ [44], though the many control points are still technically optimization parameters. Transformations are entirely described by the positions of the control points relative to their initial location before registration. Optimization of such a large number of parameters can be computationally intensive, but provides for a much larger solution search space. If the number of control points used is below the number of pixels or voxels in the image then interpolation can be used to define the transformation at more locations. The use of different types of interpolation does not tend to significantly affect results since there are many closely-spaced control points and the interpolating functions have local support.

Many groups have investigated the use of diffeomorphisms for their transformations [3, 23, 51, 28]. Diffeomorphisms are invertible transformations where both the transformation and its inverse are differentiable. Consideration of diffeomorphic transformations can help to increase registration speed and ensure that all output transformations are physical (invertible). However, this approach can become restrictive in the presence of sliding organ boundaries which can correspond to undifferentiable transformations [44].

23


1.5.4

Optimization

As previously alluded to, optimization algorithms are required to find a transformation that minimizes the cost function. Several common approaches include gradient descent, Gauss-Newton and the limited memory Broyden-Fletcher-Goldfarb-Shanno methods (lBFGS) [46, 44]. The gradient descent method finds an improved estimate of the transformation by following the gradient in the cost function. Though computationally efficient, this method is not particularly robust as it is prone to being caught in local minima. The Gauss-Newton method involves a modification of Newton’s method. Newton’s method approximates the value of cost function for the candidate solution (transformation y) plus some potential perturbation (s) by a quadratic Taylor expansion. 1 J(y + s) = J(y) + s∇J(y) + s∇2 J(y)s 2

(1.27)

where ∇2 J is like ∇J but contains the 2nd order derivatives. By setting the derivative of this expression to zero, the optimal perturbation (with respect to minimizing the cost function) can be determined. s=− (∇2 J)jk = 2

∇J ∇2 J

(1.28)

∑ ∂Ji ∂Ji i

∂yj ∂yk

∇J = 2

∑ i

Ji

+ Ji

∂Ji ∂yi

∂ 2 Ji ∂yj ∂yk

(1.29)

(1.30)

The Gauss-Newton method drops the higher order term in the expression for the components of ∇2 J. Computaton of these terms can be computationally demanding and so the Gauss-Newton method can increase optimization speed. A line search is also performed in the direction of the optimal perturbation to find a perturbation that results in a sufficiently large reduction in J. The BFGS optimization scheme uses a similar Taylor series approach but uses an

24

Chapter 1. Introduction alternative method for approximating the Hessian, ∇2 J. H = ∇J(yk+1 ) − ∇J(yk ).

(1.31)

yk corresponds to the transformation from the previous iteration of the optimizer and yk+1 is the transformation from the current iteration. The limited memory implementation of the BFGS method uses the same algorithm but involves less memory through selective storage of intermediate optimization variables. Consequently the l-BFGS method is known to be an efficient choice when working with large datasets [44]. Multi-resolution registration methods are a popular approach towards smoothing optimization and increasing speed [44, 54]. Registration is first performed on a low resolution version of the images. The transformation output by this registration is then used to initiate the registration of the images at a slightly higher resolution. This process continues in a coarse to fine resolution manner, until reaching a desired final optimization resolution. This arrangement helps the optimization procedure to avoid local minima and increases speed since each registration is initiated with the optimal transformation from the previous registration.

Figure 1.8: Multi-resolution representation of an abdominal DCE-MRI image.


1.6 1.6.1

25

Registration of Contrast-Enhanced Data Problem

A primary assumption of common registration distance measures is that changes in the images are solely caused by the motion of the patient. This assumption is no longer valid within the context of contrast-enhanced imaging. An extreme example of this is illustrated by the SSD distance measure. The SSD assumes that the pixel intensities of tissues remain constant but that they may move spatially. In the presence of contrastenhancement the intensities vary and move spatially. Consequently, the SSD distance measure does not solely reflect the degree of overlap of the images and cannot inform the registration correctly. Unlike the SSD, the MI and NMI distance measures do not assume that pixel intensities remain constant between images. However, these distance measures are also affected by contrast-enhancement. Successful registration is dependent on increasing the value of the joint probability density function (jPDF). This corresponds to improving the ability to predict the intensities in one image, if the intensities in the other image are known. MI based registration relies on an assumption of the existence of a one-to-one mapping of the intensities in one image to the intensities in the other. The problem lies in the partial enhancement of structures within the images. Consider an image of an unenhanced organ that has a uniform intensity of 1. At some later time, the organ has not moved, but half of it has enhanced to an intensity of 2 while the other half has enhanced to an intensity of 3 (Figure 1.9). Although the two images are perfectly aligned, it is not possible to predict with certainty the intensities in the enhanced image if the intensities in the unenhanced image are known. A one-to-one mapping of the intensities does not exist which results in a value of the jPDF, and consequently the distance measure, that is unrelated to motion. The NCC distance measure allows for intensity differences between the images to be


26

Figure 1.9: Mock image of a kidney-like organ with partial enhancement after time ∆t. registered but assumes linear dependence. Since the location and magnitude of contrast enhancement is dependent on local vasculature and flow rates, enhancement patterns tend to be very heterogeneous within contrast-enhanced images. Structures in the images can partially enhance or enhance at different times therefore the enhancement patterns cannot be adequately described by linear relationships. Consequently, the NCC distance measure is not well suited to registration of contrast-enhanced data. The NGF distance measure is also affected by contrast-enhancement. NGF assumes that the same gradients occur within a pair of images despite the images having different intensities or coming from different modalities. The presence of the contrast agent can significantly increase image detail and may produce new image gradients that are not present in the unenhanced image. Successful registration relies on finding a transformation that matches the locations of the gradients in the two images. Since many of the gradient locations cannot be matched within this scenario, the quality of registration may be reduced. Ultimately, changes related to contrast enhancement must somehow be differentiated from changes due to motion. In addition, patient motion must be removed while simultaneously preserving the contrast-agent dynamics that are the focus of subsequent analyses. Consequently, motion correction for contrast-enhanced datasets presents a challenging problem that remains a very active area of research [12, 39, 42, 68, 1, 18, 19, 21, 27, 37, 40, 41, 52, 57, 62, 67].


1.6.2

27

Approaches

Relatively few registration techniques have been developed within the context of abdominal and internal thoracic DCE-MRI [25, 12, 42, 39, 5, 68, 62]. These sites pose a challenging registration problem due to the large amounts of deformable motion induced by respiration. Several approaches have involved registering to a set of constructed static reference images that contain time-varying contrast enhancement [12, 42, 62]. Each image in the dataset is registered to its own reference which has been constructed such that it contains a similar level of contrast enhancement. However, these approaches tend to rely upon having a temporally well-sampled dataset and therefore may not be well suited to registration of sparse data (e.g. pair wise registration). Buonaccorsi et al [12] used pharmacokinetic modelling to simulate time-varying enhancement within a static reference volume. The model is initially fit to the unregistered data. After registration, new reference volumes are constructed by fitting the pharmacokinetic model to the registered data. An iterative process of reference generation and registration occurs until stopping criteria are met. Registration was simplified by considering only translational motion, which corresponds to an assumption of tissue rigidity. Since tumours tend to be more rigid than normal tissue, the distance measure was calculated solely using the tumour ROIs. Melbourne et al [42] proposed a nonrigid technique that does not rely upon the use of models or ROI definitions. Principal components analysis (PCA) reconstruction (Appendix B) was used to create a set of enhancing reference images from the unregistered dataset. Similar to Buonaccorsi et al [12], a new set of reference images were created using the newly registered data, followed by an iterative process of reference generation and registration. Additional components are included in the PCA reconstruction as registration continues in order to incorporate more complex uptake profiles into the reference image construction. PCA based motion correction relies upon the assumption that motion is largely uncorrelated in time. Uncorrelated motion is removed by reconstruct-


28

ing using only the principal components which contain the most temporally correlated changes in the data. Consequently this technique was observed to fail when periodic motion appeared in the principal components used for reconstruction [42]. Wollny et al [62] developed a unique registration algorithm that exploits the quasiperiodicity within free-breathing myocardial MRI perfusion scans. Similar to [42], the technique is nonrigid and does not use pharamacokinetic models or ROI definitions. Sequence similarity profiles are used to select a subset of image volumes that correspond to a similar phase of the respiratory cycle. The subset is registered together using the NGF distance measure and then treated as as a static yet temporally down sampled version of the original perfusion scan. A reference is generated for a given image within the original DCE-MRI dataset by interpolating the temporally down sampled scan at the time of acquisition of the given image. The algorithm was demonstrated to successfully register clinical MRI perfusion scans and remove motion artifacts from intensity-time curves. However, since registration was only performed in two dimensions, reduced accuracy was observed in the presence of out-of-plane motion. Furthermore, difficulties were encountered producing a temporally well-distributed subset when a patient breathed in a shallow manner. This was observed to reduce registration accuracy. The assumption of quasi-periodicity may not be broadly applicable for other sites or within imaging contexts with lower temporal resolutions.


1.7

29

Hypothesis

DCE-MRI is an important tool for the investigation of the vascular environment within many different sites [45, 66, 47, 68, 42, 53, 24, 2, 50, 49, 34, 15, 6]. Involuntary patient motion (e.g. respiration) is unavoidable within many of these imaging contexts and can significantly decrease the accuracy of analysis to the point where the data can be rendered unusable [48]. Image registration has effectively removed motion within many imaging contexts but is challenging to develop for contrast-enhanced imaging. Consequently, registration for contrast-enhanced data remains an open problem. Current approaches often take advantage of unique characteristics of their imaging site and make generally reasonable assumptions about the image data. Common assumptions include rigid, non-periodic, or periodic motion, as well as tissue incompressibility [12, 42, 62, 68]. Though favourable registration results are observed in many cases, poor results usually relate to special cases where initial assumptions are violated. The ideal registration technique should make minimal assumptions about the image data in order to be robust and widely applicable. This investigation presents and tests a novel fully-automated approach for the registration of DCE-MRI data. It is hypothesized that an intensity correction derived directly from the image data can facilitate registration without the use of assumptions related to motion characteristics or the use of pharmacokinetic models.

Chapter 2 Registration Using Motion-Informed Intensity Corrections The proceeding investigations have been submitted in modified form to Magnetic Resonance in Medicine and are described in further detail here.

30

Chapter 2. Registration Using Intensity Corrections

2.1

31

Introduction

Here we propose an alternative registration approach that uses a novel motion-informed intensity correction to facilitate the registration of DCE-MRI datasets. The intensity correction simulates the presence or absence of the contrast agent in the image volumes to be registered in an attempt to emulate the level of contrast-enhancement present in a single reference image volume. Intensity differences due to the contrast agent are reduced resulting in a reassertion of the assumption of patient motion as the sole source of changes in the images to be registered. The resulting deformations are then applied to the original data. This nonrigid 3D registration approach does not require the use of pharmacokinetic modelling or ROI definition and should be robust in the presence of periodic motion. In addition, the technique is unrelated to the temporal sampling of the image data. The algorithm is applied to clinical DCE-MRI datasets and evaluated quantitatively using a simulation test case. The potential clinical utility of the registration technique is investigated through analysis of pharmacokinetic modelling performed before and after registration.

2.2 2.2.1

Materials and Methods Patient Data

3D abdominal DCE-MRI scans of 10 patients with advanced renal cell carcinoma were utilized for this study. The patients were scanned as part of an ongoing clinical trial at our institution to investigate the capability of DCE-MRI, DCE-CT, and DCE-US to assess tumour response to the anti-angiogenic drug sunitinib malate (Sutent). 3D scans were acquired using a fast spoiled gradient-echo sequence on a 1.5 Tesla General Electric scanner. The spatial resolution was 1.88 mm by 1.88 mm by 8 mm in the LR, SI,


32

and AP directions respectively, with each image volume consisting of 256 by 256 by 12 voxels. The temporal resolution was approximately 3.7 seconds with a total of 80 image volumes acquired for each patient scan. The contrast agent gadodiamide (Omniscan, GE Healthcare) was injected 20 s into the scan. The patients were scanned before and after treatment. A radiologist delineated the boundaries of the tumour on a single image volume within each scan.

2.2.2

Image Registration Framework

Image registration was framed in terms of an optimization problem. Consider the task of finding a nonrigid transformation y that registers some moved image volume T to a reference image volume R. The transformation is found by minimizing the objective function J(y) = D(T (y), R) + αS(y)

(2.1)

where T (y) corresponds to T transformed by y. D is a distance measure which quantifies the similarity between T (y) and R, while S(y) is a regularizing term which quantifies the suitability of the candidate transformation y. α is a scalar weighting parameter to control the degree of regularization. Here the distance measure has been chosen to be the common sum of squared differences of the voxel intensities (SSD). Traditionally the SSD performs poorly in the presence of contrast-enhancement due to its high sensitivity to intensity differences. However, it is precisely this sensitivity that is taken advantage of by the intensity correction utilized within this work. The regularizer was chosen to be the elastic regularizer [11], with [µ, λ] = [1, 0] and α = 800. A non-parametric transformation model is employed [44]. Transformations are described by a grid of control points that are superimposed on the image volume to be registered. Displacement of the control points displaces the corresponding regions of the image. Registration is achieved by finding a set of control point positions which describe


33

a transformation that aligns the features of T with the features of R. Multi-resolution optimization is used in order to smooth the optimization problem and increase speed [44]. An initial rigid registration is performed prior to multi-resolution optimization to correct for global motion. The optimizer was selected to be the limited memory Broyden-Fletcher-Goldfarb-Shanno method which is well suited to working with large datasets [44]. The Flexible Algorithms for Image Registration (FAIR) toolkit [44] was used to implement this framework.

2.2.3

Reference Selection

The reference image volume R was chosen to be the volume with the maximum average voxel intensity. This volume is similar to the majority of the other image volumes within the dataset in terms of the amount of contrast enhancement. This presents an easier registration problem than the more traditional method of choosing the reference to be a pre-enhancement volume which has significantly different image intensities. In addition, the contrast agent provides additional image detail which provides better local information for the registration algorithm to use.

2.2.4

Intensity Correction

A primary assumption of intensity-based registration using the SSD distance measure is that intensity changes in the image volumes are only due to motion.

However,

with contrast-enhanced datasets this assumption is no longer valid since image intensity changes can result from both motion and the presence of the contrast agent. In an attempt to restore the validity of this assumption an intensity correction (c) is applied to the moved image volume T to account for differences between T and R that are caused by the presence or absence of contrast agent. The intensity corrected volume (Tc ) and reference image volume R are then registered together and the output transformation is applied to the uncorrected image volume T . The intensity correction term serves to re-


34

duce the contribution of intensity differences to the SSD distance measure that are due to the contrast agent. As a result, un-enhancing regions of the image play an increased role in informing the image registration. This is desirable since registration of such regions is a simpler task. Ideally, this correction should alter the appearance of the structures in T to have the same amount of contrast enhancement as R. The difference between R and T serves as a starting point for finding such a correction. Given that these two image volumes were perfectly aligned, R − T would be precisely the correction we are looking for. When the structures in the two image volumes are not perfectly aligned, R − T becomes an approximation to the ideal correction we need to apply to T . The upper bound on the spatial uncertainty of R − T as an approximation to the ideal correction is given by the maximum amount of motion present between R and T . However, R − T cannot be directly used as an intensity correction since,

Tc = T + c = T + (R − T ) = R

(2.2)

The problem relates to the effective spatial uncertainty of R − T when used directly for correction. R−T contains information specified at the resolution of the images and defines a unique value for the intensity correction at each pixel (voxel). However, the spatial uncertainty of R − T as an intensity correction is equal to the magnitude of motion, which exceeds the size of a pixel (voxel). Therefore if it is applied unaltered we are assuming higher spatial precision than is actually supported by the inherent uncertainty. A somewhat analogous example would be the reporting of a number with more precision than the number of significant digits. Therefore the uncertainty is not explicitly present, rather it is only implied by the way in which we seek to use R − T as an intensity correction. In order to ensure that the intensity correction is in effect only as spatially accurate as this uncertainty limit, we

35


Figure 2.1: Intensity correction for a simulated enhancing necrotic tumour. blur R − T using a Gaussian blurring kernel. We define the intensity correction to be Tc = T + c

(2.3)

c(y) = (R − T (y)) ∗ N (0, σ)

(2.4)

where N (0, σ) is a 3-dimensional Gaussian blurring kernel with zero mean and standard deviation σ. The uncertainty limit can then be imposed by setting the full-width half maximum (FWHM) of the blurring kernel to be equal to twice the maximum magnitude of motion present between R and T . The factor of two is included to allow for motion in both directions along each coordinate axis. With each subsequent iteration of the registration algorithm a new and potentially improved candidate transformation y is considered for alignment of T and R. The intensity correction is recalculated using T transformed by the current candidate transformation (Equation 2.4). Therefore the intensity correction is recalculated upon each iteration of the registration algorithm. The same blurring kernel is used every time the intensity correction is recalculated. Altering the blurring kernel size during registration is the subject of further discussion in section 3.2.1. An intensity correction example is illustrated in figure 2.1. Here, a simulated necrotic tumour moves and enhances. The unenhanced image serves as the moved image T , and the enhanced image is the reference R. The difference between the two images is blurred (|c|) and added to T to form Tc . The absolute difference between T and R serves as a reasonable surrogate for investigating how the intensity correction affects the SSD. In figure 2.2, T , R, and |R − T |


36

Figure 2.2: Simulated tumour images I) before intensity correction II) after intensity correction III) after an ideal intensity correction

are illustrated before intensity correction (I), after intensity correction (II), and after an artificially constructed ideal intensity correction (III). The ideal correction eliminates all of the differences between the images due to enhancement. In III, the SSD surrogate is only dependent on the motion of the tumour whereas in I the surrogate is dependent on both motion and enhancement. Upon intensity correction (II), the effects of enhancement on the surrogate are reduced, and it takes on characteristics that are similar to III. Accurate estimation of the maximum amount of motion present between an arbitrary pair of image volumes within a contrast-enhanced dataset poses a challenging problem. For practical implementation purposes this must be an automated process. We use an alternative approach which involves determining the maximum amount of motion present within the entire DCE-MRI scan (Mmax ), rather than individually for each pair of image volumes within the scan. This single global maximum is then used to define the FWHM of the blurring kernel used during the registration of each pair of image volumes within


37

the scan. LR SI AP F W HM = 2[Mmax , Mmax , Mmax ]

2.2.5

(2.5)

Estimating Mmax

Mmax is estimated for a DCE-MRI scan by taking advantage of the slower intensity changes that occur during contrast agent wash-out. Registration is much easier in the wash-out phase than in the wash-in since intensity changes are much less drastic. Although conventional registration techniques are generally still insufficient for wash-out volumes, they are potentially useful for estimation of Mmax . To estimate Mmax , 5 maximally dissimilar pairs of image volumes are identified from the wash-out phase using the SSD distance measure. A temporal windowing technique is used while searching for the dissimilar pairs in order to ensure that they are dissimilar because of motion rather than contrast agent wash-out. For example, a maximally dissimilar pair comes from a window that includes a maximum of 10 volumes acquired consecutively in time. The aforementioned registration technique is then used to register the pairs to one another, using increased regularization (α = 1000) and without intensity correction. The distributions of voxel displacements for each of the 5 registrations are analyzed and the 99th quantile is selected to represent the maximum amount of motion for each registration. The maximum of these 5 estimates corresponds to Mmax .

2.3

Evaluation

The evaluation of the proposed registration algorithm involves two main elements. First, the accuracy of the registration is assessed quantitatively using a DCE-MRI dataset with simulated motion. Second, potential clinical impact is investigated by studying the effect of image registration on the fitting of a pharamacokinetic model to ROI intensity-time curves from multiple clinical DCE-MRI datasets.


2.3.1

38

Simulation

Ideally, validation would involve simulating a known amount of motion on a motionless clinical DCE-MRI dataset. Image registration would then be applied and the distance of voxels from their initial location in the motionless dataset would be compared before and after registration. In practice, there is always some motion present in DCE-MRI scans so this is not possible. Consequently the motionless dataset must be constructed. One approach could involve simulating time-varying contrast-enhancement on a single pre-contrast image volume. Though this method ensures that the dataset is completely motionless, simulated contrast-enhancement may not accurately reflect the complex characteristics of real enhancement patterns. These characteristics make DCE-MRI registration a challenging problem so we have opted to preserve them within the motionless dataset used for validation. To produce a motionless dataset, a two step approach has been taken. First, the aforementioned image registration technique is applied to a clinical dataset to remove a large amount of the motion present. Second, a principal components analysis reconstruction of the registered data is performed to remove remaining temporally-uncorrelated motion (Appendix B). It is assumed that the majority of the temporally correlated motion (e.g. respiration related) is removed or becomes uncorrelated due to the initial registration. The first four principal components are smoothed and then used for reconstruction. After constructing the motionless dataset, realistic deformable motion must be simulated. When a patient inhales while lying supine, structures in the abdomen tend to move in the anterior and inferior directions [8, 13]. In an attempt to simulate this motion, the diaphragm has been modeled as a half sinusoid with temporally varying amplitude. Let x, y, and z correspond to the LR, AP, and SI coordinate directions respectively. Consider a point P located at an arbitrary position (x0 , y0 , z0 ) within the motionless dataset. Within the dataset with simulated motion, the position of P is dependent on x0 and the time t that each image volume was acquired. The position of P within the simulation is


39

described by, x(x0 , t) = x0

(2.6)

y(x0 , t) = y0 + ∆y(x0 , t)

(2.7)

z(x0 , t) = z0 + ∆z(x0 , t)

(2.8)

∆y(x0 , t) = ∆APmax |sin(πt/tb )|sin(πx0 /xmax )

(2.9)

∆z(x0 , t) = −∆SImax |sin(πt/tb )|sin(πx0 /xmax )

(2.10)

where

∆SImax and ∆APmax correspond to the maximum SI and AP displacements respectively. xmax is the maximum LR extent of the patient, and the respiration rate is given by 1/tb . The respiration rate is sampled from a Gaussian distribution upon each new breath and is centered about the average respiratory rate of humans, 15 breaths/min (σ = 2 breaths/min) [4]. Finally, a small rotation about the z (SI) axis is added, where the angle of rotation is directly related to the magnitude of the simulated breath. The angle of rotation at time t is given by, θ(t) = |sin(πt/tb )|θmax

(2.11)

where θmax = 1 °is the maximum permitted rotation. The rotation introduces some additional displacement of structures in LR and AP directions. These transformations are applied to the motionless DCE-MRI dataset to obtain a DCE-MRI dataset with simulated motion. This dataset is subsequently registered using the proposed intensity correction based registration technique. The transformation output by this registration is compared to the known inverse of the simulated motion to establish quantitative accuracy. The accuracy of competing registration techniques can also be determined through registration of this dataset. The accuracy of the proposed and competing techniques can then be compared to help illustrate efficacy. The full validation pipeline is summarized in Figure 2.3.


Figure 2.3: Quantitative validation pipeline

40


41

A control dataset is also constructed in order to estimate an upper bound for potential registration accuracy and lend context to the quality of registration results. The control is created by simulating an identical sequence of motion in the first pre-enhancement volume of the motionless dataset. This produces a dataset with motion characteristics that are identical to those of the previously described simulation (equations 2.10, 2.9, and 2.11) but without contrast-enhancement. Since the proposed intensity correction seeks to eliminate the effects of contrast agent related intensity changes, the maximum registration accuracy should be well estimated by registering the un-enhancing control dataset. Registration of the control dataset is performed using the the SSD, elastic regularizer, and multi-level optimization but without using intensity correction.

2.3.2


A modified 2-compartment Tofts model is used to determine ve and K trans parameters for the renal tumours within the abdominal DCE-MRI data [59]. The arterial input function is determined using the signal from a small region of descending aorta near the renal tumour. The concentration of the contrast-agent in the blood plasma is determined for this region at each time point in the scan. This produces an estimate of concentration that is explicitly dependent on the image data at all time points. This is in contrast to the unmodified Tofts model where only the initial value of the concentration is dependent on the image data. With the unmodified model, the AIF is assumed to be a delta function with a height determined from the image data at a single initial point in time. All subsequent dynamics are then determined by the model. The assumption of a delta function AIF is a useful simplification, however it may not accurately represent the dynamics in the blood plasma. Measuring the concentration of the contrast agent in the blood plasma at all times ensures a more accurate representation of these dynamics. The modified model is fit to the pre and post-treatment scans of the 10 patients, before and after registration, using the average signal from the tumour ROIs delineated


42

by the radiologist. In general, the tumour ROI for each patient was delineated on a single image volume that contained maximum contrast enhancement. Mean value ROI analysis is used as it is the analysis protocol adopted within the clinical trial that is underway at our institution. A coarse voxel-based model fitting is also performed for a single dataset in order to investigate the effects of registration on a more local analysis of the tumour. The in-plane resolution of the dataset is down-sampled from 1.88 mm by 1.88 mm to 7.52 mm by 7.52 mm. Model fitting is then performed independently for each voxel producing local ve and K trans values. The resulting parameter maps are then linearly up-sampled to the original 1.88 mm by 1.88 mm in plane resolution for comparative purposes.


2.4 2.4.1

43

Results Simulation

A DCE-MRI dataset with complex enhancement patterns and no visually assessable motion was successfully constructed using the method outlined in the evaluation section. Motion was then simulated using ∆SImax = 19.7mm and ∆APmax = 6.8mm. These two values correspond to the average SI and AP components of the 20 Mmax estimations that were determined during registration of the clinical DCE-MRI datasets. The dataset containing realistic contrast enhancement and simulated motion was registered using the proposed method and an alternative nonrigid registration technique that employs the NGF distance measure combined with multi-level optimization and elastic regularization [44]. The NGF distance measure establishes correspondences between images volumes based on the location of image gradients (and not gradient strength) and is therefore expected to be more robust than conventional registration techniques in the presence of intensity changes. The proposed registration algorithm reduced the average motion in the simulated motion dataset from 8.9 mm (σ = 3.6mm) to 2.9 mm (σ = 1.1mm). The NGF-based algorithm was less successful, reducing the average motion to 4.9 mm (σ = 1.8mm). As expected, the control dataset was most accurately registered with 2.0 mm (σ = 0.8mm) of motion remaining after registration. Figure 2.4 indicates the distribution of the distances of voxels from their original locations before and after these registrations. For clarity of presentation, only the registration accuracies associated with the first 40 image volumes are shown. These image volumes represent the most challenging registration problems within the dataset as they contain the most significant intensity changes. Therefore the accuracy of the proposed registration algorithm is well-demonstrated by these registration accuracies. A modest decrease in the accuracy of the proposed method was observed for pre-enhancement and contrast agent wash-in volumes (Figure 2.4c, time < 74s).


44

Figure 2.4: Distribution of the distance of voxels from a reference image volume before registration (a), after registration using NGF (b), after registration using the proposed method (c), and after registration of the control (d) for the first 40 time points. Box plots encompass the 25th to 75th quartiles and whiskers encompass 97.5% of all values. The average is denoted by the black line.


45

Figure 2.5: Relative changes in chi-squared scores for the Toft’s model fits due to registration. Scores were normalized by the pre-registration values to indicate relative changes for each scan. The average change (AVG) in the relative values is shown for the set of chi-squared scores that increased and decreased. Values under unity correspond to an improved fit.

2.4.2


The modified 2-compartment Toft’s model was fit to the pre and post-treatment DCEMRI scans, before and after registration. The reduced chi-squared scores were calculated in order to determine whether registration improved model fitting (Figure 2.5). The chi-squared scores decreased after registration for 7/10 cases within both the pre and post-treatment sets of scans, indicating improved model fits. A small increase in the chi-squared scores was present for 3/10 cases within both the pre and post treatment sets of scans indicating slightly worse fits. Two of the three unimproved fits involved the same two patients for both the pre and post-treatment scans. Upon visual assessment, limited tumour motion was found within the scans of these two patients before registration was performed. The other unimproved scans contained more homogeneous tumours and surrounding tissue that exhibited similar enhancement patterns. Consequently, the effect of removing motion did not significantly alter the average signal used for model fitting leading to a slightly worse model fit.


46

Figure 2.6: Distribution of chi-squared scores after registration for the coarse voxelbased fits (normalized by pre-registration scores). Box plots encompass the 25th to 75th quartiles. Whiskers encompass 97.5% of all values. Values under unity correspond to an improved fit. Voxel-based fitting was performed for a single dataset before and after registration. The four ROIs that were delineated by the radiologist contained an average of 264 voxels each. The reduced chi-squared scores for each of the voxel fits were calculated before and after registration. Similar to the mean ROI analysis cases, significant reductions in the normalized chi-squared scores of the voxel-based fits were observed (Figure 2.6). Parameter maps were generated for the four ROIs in order to investigate the local changes in parameter estimation due to registration (Figures 2.7 and 2.8). Motionrelated estimation artifacts reduced for both parameters after registration. Improvements in spatial localization are observable in the K trans map after registration, predominantly near regions with high K trans values. Prior to registration, these regions of tissue would move back and forth, contaminating the signal of adjacent voxels. Consequently adjacent voxels take on characteristics that are associated with higher K trans values that may not actually reflect their true value. This produces a small washing out effect that manifests itself in the pre-registration maps. For the ve map, the most notable improvements are in the anterior portions of the tumour, next to the tumour-liver boundary. Prior to


47

Figure 2.7: K trans parameter maps for ROIs from a single pre-treatment scan, before (top) and after registration (bottom). Peak enhancement images of the ROIs are shown in grayscale. Arrows indicate major regions of improvement. Mean-value ROI analysis K trans values for ROI 1-4 are [0.085, 0.12, 0.16, 0.50] min−1 and [0.070, 0.084, 0.10, 0.21] min−1 before and after registration respectively. registration, the highly perfused segment of liver moves in and out of the regions being investigated which contributes to an overestimation of ve . After registration, the signal contamination caused by the liver is eliminated resulting in more suitable ve estimates.


48

Figure 2.8: ve parameter maps for ROIs from a single pre-treatment scan, before (top) and after registration (bottom). Peak enhancement images of the ROIs are shown in grayscale. Arrows indicate major regions of improvement. Mean-value ROI analysis ve values for ROI 1-4 are [0.13, 0.14, 0.21, 0.25] and [0.085, 0.09, 0.11, 0.14] before and after registration respectively.

49


350 before registration after registration

Voxel Intensity

300 250 200 150 100 50 0

75

150

225

300

Time(s)

Figure 2.9: Intensity-time curve before and after registration for a voxel with an improved pharmacokinetic model fit after registration. A single voxel located within the hyperintense region of the ve parameter map for ROI 1 in Figure 2.8 (indicated by a white arrow) was selected to demonstrate how registration improves intensity-time curves. Prior to registration the intensity-time curve for the voxel is severely contaminated by signal from adjacent tissue and consequently does not exhibit the typical phases of contrast-enhancement: pre-injection, wash-in and wash-out (Figure 2.9). This reduces the quality of the pharmacokinetic model fit and results in an unphysical inflated local estimate of ve . After registration, signal contamination is reduced and the intensity-time curve displays more typical enhancement characteristics resulting in an improved model fit and value of ve .


2.5

50

Discussion

The registration algorithm successfully removed the majority of the visually assessable motion within all 20 clinical datasets. Estimations of Mmax appeared to correlate well with the magnitude of motion present in different scans. The average LR, SI and AP components of the 20 Mmax estimations used for the registration of the clinical DCEMRI datasets were 7 mm, 20 mm, and 7 mm respectively. These estimates agree with typical values for abdominal organ motion found in literature [8, 13]. The estimation of Mmax for the simulation was 3 mm, 21 mm, and 7 mm in the LR, SI, and AP directions respectively. These values compare well with the true values of 0.2 mm, 19.7 mm, and 6.8 mm indicating that the Mmax estimation algorithm produced reasonable results. Average motion within the simulation was shown to be significantly reduced (2.9 mm, σ = 1.1mm) and compared well with theoretical ’best-case’ accuracy limit implied by registration of the control (2.0 mm, σ = 0.8mm). The proposed method also out performed a competing registration algorithm based on the NGF distance measure. Ultimately, an ideal registration algorithm endeavours to align contrast-enhanced data in a way that is completely independent of the contrast-enhancement. The control dataset has no enhancement and has anatomical and motion characteristics that are identical to the dataset that does enhance. Consequently, registration of the control lends insight into how such an ideal algorithm might perform. Registration of unenhancing data is a much simpler task for which many successful registration techniques exist. As a result, the associated registration accuracy sets a very high standard that should be greater than or equal to the accuracy of competing contrast-enhanced registration techniques. Stateof-the-art techniques are often complex, site specific, and difficult to reproduce properly without direct collaboration with their developers. Consequently comparisons with such techniques appear infrequently in the literature. In this work we have compared our method with a competitive algorithm based on the NGF distance measure. We selected this method because the NGF is expected to perform well in the presence of contrast


51

enhancement and a 3D implementation was available to us [44]. However, registration of an unenhancing control dataset using conventional techniques may provide a more accessible comparative alternative and set a greater standard of accuracy. Consideration of a finer resolution by the multi-level optimization scheme should improve accuracy at the cost of computational speed. The highest resolution considered within this study was 7.52 mm by 7.52 mm by 3mm in the LR, SI, and AP directions respectively. The actual resolution of the datasets in this study is 1.88 mm by 1.88 mm by 8 mm. Consideration of a higher final resolution would allow the algorithm access to more detailed local image information within the coronal plane which could be used to make finer local adjustments to the deformation grids. A graphics processing unit implementation of the FAIR toolkit would facilitate the consideration of finer resolution levels by increasing computation speed. The uniform use of a single maximally-sized blurring kernel for all pair wise registrations within a dataset provided for robust registration results. In principle, pairs of image volumes with less motion than Mmax could be registered using a smaller blurring kernel. Since favourable results were still achieved for volumes with significantly less motion than Mmax , the proposed registration technique appears to be insensitive to overestimation of the size of the blurring kernel. Therefore successful registration may often be achievable using ”best guess” estimates of Mmax as long as the true value is exceeded. This ensures the utility of the algorithm in contexts where Mmax cannot be estimated using the previously described method. Despite these findings, the use of an unnecessarily large blurring kernel could potentially have an adverse effect on registration quality. When the blurring kernel is too large, the spatial uncertainty in the intensity correction exceeds its true value which could degrade registration results. Ideally a uniquely sized blurring kernel would be used for the registration of each pair of image volumes. In this way the spatial uncertainty of the intensity correction would be as low as is achievable, potentially improving registration


52

results. Further gains could be realized if the maximum motion between a pair of image volumes could be estimated dynamically during registration. The blurring kernel size would reduce as registration improved the alignment of the volumes. This would result in an improved estimate of the intensity correction, which would in turn further facilitate registration. Dynamic estimation of the maximum motion present between a single pair of image volumes will be the focus of future investigations. In order to create the simulation test case, a motionless dataset needed to be constructed from a clinical DCE-MRI scan. The proposed registration technique was used first and then motion was further eliminated using PCA reconstruction methods. Whether or not this combined method for motion correction is valid outside of a simulation context is an open question. Such a method would be difficult to validate since PCA reconstruction does not output a deformation which can be compared to some ground truth. Furthermore, PCA reconstruction necessarily omits some of the image content in order to eliminate motion. This could have an adverse effect on subsequent DCE-MRI analysis if some of the contrast agent dynamics were unknowingly omitted during reconstruction. This approach may have utility however by providing a set of static reference images with time-varying enhancement. The source DCE-MRI dataset could then be registered to the reference images, similar to the methods proposed by Buonaccorsi and Melbourne [12, 42]. Image registration was shown to improve mean-value pharmacokinetic model fits as indicated by decreases in the reduced chi-squared scores. Although some fits were unimproved, these generally corresponded to scans with either highly homogeneous tumours or limited tumour motion prior to registration. These results suggest that image registration has an important role to play in improving mean value ROI analysis even though it is more robust to motion than voxel-based analysis. The improvement in voxel-based DCE-MRI model fits was naturally much more apparent after image registration. The magnitude of the motion within the dataset before


53

registration far exceeded the size of the investigated voxels. Therefore contamination of signal from tissue in adjacent voxels occurs before registration as the underlying tissue moves. Registration successfully reduced fitting artifacts related to this issue. However, as the underlying tumour becomes more homogeneously enhancing, the negative effects of motion on voxel-based ve estimates become less observable. This is indicated in figure 2.8 where we see fewer observable fitting artifacts for ROIs 3 and 4 before registration. Despite this, local trends become clearer for these ROIs after registration as is indicated by improved spatial correspondence between ve values and the peak enhancement images. Mean value ROI analysis was found to be insensitive to small regions within the heterogeneous tumour investigated using voxel-based analysis. Average post-registration K trans and ve values were between 0.07-0.21 min−1 and 0.085-0.14 respectively. These are much lower than the non-zero values present in the parameter maps due to signal averaging that takes place over the tumour’s necrotic core. Consequently, treatmentinduced changes to the tumour vasculature may be easier to identify using a voxel-based analysis of registered DCE-MRI data.


2.6

54

Conclusions

A nonrigid DCE-MRI registration technique that uses a novel motion-informed intensity correction was developed and applied to 20 abdominal clinical datasets and a preliminary validation test case. Significant motion correction was observed visually for clinical datasets and quantitatively for the validation test case. The proposed method was also shown to be more accurate than a competing method. The potential clinical utility of the registration technique was demonstrated by the post-registration improvement of pharmacokinetic modelling for both mean value and voxel-based DCE-MRI analyses. The intensity correction based registration technique should be easy to implement using other registration frameworks since it only relies upon altering the underlying image data used during registration. Dynamic motion estimation and algorithm implementation on a graphics processing unit are predicted to improve registration performance and will be the subject of future investigations.

Chapter 3 Conclusions and Future Work

55

Chapter 3. Conclusions and Future Work

3.1

56

Significance

The problems associated with registration of contrast-enhanced image data using conventional intensity-based registration techniques have been discussed. A novel registration method has been developed and shown to effectively remove motion from data with simulated motion and improve the fitting of a pharmacokinetic model to clinical DCE-MRI data. Consequently the technique may be useful for improving the quality of DCE-MRI analysis. Implementation of the method within other registration frameworks should be straightforward as it only relies upon a small manipulation of the image data during registration. Therefore the method is highly accessible and may be amenable to combination with other registration techniques to achieve improved results. Since the method does not make use of characteristics of the anatomical site or the imaging modality, it may be useful for registration of contrast-enhanced data acquired using other modalities such as DCE-CT and DCE-US. The quantitative validation approach employed within this work may be useful for lending context to future studies regarding registration of contrast-enhanced data. It is difficult to properly reproduce and optimize competing algorithms for comparison and is therefore rarely done in the literature. If quantitative validation is performed using contrast-enhanced data with simulated motion, then comparison with registration of an unehancing control dataset with identical motion and anatomy offers a compelling alternative. Registration of unenhancing data is a much simpler task and permits the use of conventional registration algorithms that have been highly refined, packaged, and are widely distributed. An ideal registration algorithm for contrast-enhanced data endeavours to be unaffected by the contrast-enhancement. Therefore comparison with registration of the unenhancing control presents a sensible alternative to reproducing and applying complex competing registration algorithms that may be highly sensitive to parameters, anatomical site, and imaging modality.


3.2

57

Future Work

This investigation has endeavoured to develop and test a new method for registration of contrast-enhanced data. Although utility has been demonstrated using both simulated and clinical examples, there are several supplementary research directions that could be explored in the future. Improved blurring kernel definition and motion estimation, increasing registration accuracy, further application and testing, as well as investigating the effects of motion on DCE-MRI analysis will be discussed.

3.2.1

Blurring Kernel Definition and Motion Estimation

A blurring kernel is used to enforce the uncertainty limit on the intensity correction that is used to facilitate registration. As previously discussed, both over and under estimation of the blurring kernel size can potentially decrease registration accuracy. In the implemented scheme, a single global maximum value for the maximum amount of motion present in the scan is estimated and then used to define the blurring kernel size for each pairwise registration. This approach ensures that the blurring kernel size is always overestimated. This method produced favourable results because registration quality was less sensitive to overestimation than underestimation of the blurring kernel size. It is hypothesized that registration accuracy could be further improved by defining a unique blurring kernel for each pair of images to be registered. This requires the estimation of the maximum amount of motion present between any two arbitrary images. In the current implementation, the global maximum amount of motion is estimated using a pre-registration of image volumes with similar levels of enhancement. This permits the use of conventional registration techniques. Since an arbitrary pair of images could differ greatly in terms of their level of enhancement, conventional registration techniques will be unable to estimate the maximum motion in some cases. Several approaches could be investigated to overcome this challenge. The registration


58

algorithm presented in this work could be used to perform a coarse pre-registration of each image volume pair. The blurring kernel size would be defined by the global maximum motion value or alternatively an educated guess that overestimates the global maximum motion. The transformations resulting from the pre-registration would then be used to define unique blurring kernel sizes for each subsequent pairwise registration. This approach would significantly increase the total runtime for the registration algorithm. Integration of the coarse pre-registration step with the subsequent full registration could eliminate this issue. Registration would proceed exactly as is detailed in Chapter 2 except the blurring kernel size would change to a unique value for each pair of images after the registration reaches a specified resolution level. This resolution would be equivalent to the level previously used for the coarse pre-registration step. Theoretically as registration proceeds, the motion present between the two image volumes decreases. Ideally, the size of the blurring kernel would decrease along with the decreases in motion. In this way, the intensity correction would become more and more spatially accurate which could further improve subsequent registration. Implementing such a scheme poses a challenging problem that involves dynamically estimating the maximum amount of motion remaining between a pair of images throughout registration. This would likely involve sub-registration steps that would considerably increase computation costs. Implementation of such a scheme on graphics processing units could decrease computation time. Changing the size of the blurring kernel mid-registration once or dynamically throughout registration could have an adverse effect on convergence. The multi-resolution scheme used in this work helps the optimizer to find a transformation that corresponds to a more global minimum in the objective function. Changing the blurring kernel size midregistration may displace the current solution to a less optimal local minimum. Such a change may be irreversible because the multi-resolution scheme becomes less likely to permit escape from local minima as it proceeds. Such effects would be avoided if the


59

registration was reset to the coarser resolution level every time the blurring kernel was redefined, however this may unacceptably increase computation time. Restrictions on the number of registration resets could be present a possible solution. Definition of the blurring kernel size is integral to the success of the proposed registration algorithm. Future algorithm improvements should focus on individualized and dynamically defined blurring kernel sizes for each pairwise registration. These avenues are liable to require graphics processing units to ensure practical computation times but may provide for increased registration accuracy.

3.2.2

Resolution and Registration Accuracy

The abdominal DCE-MRI data registered within this study had a voxel size of 1.88 mm by 1.88 mm by 8 mm in the SI, LR, and AP directions respectively. For the abdomen, the second largest component of motion is typically in the AP direction [13, 8]. The image data is sparsely represented in this direction which could potentially have an adverse effect on AP motion correction. Sparse representation could cause the appearance of local structural information in one image volume that is not present in another due to movement in the AP direction. It also results in a reduced number of unique AP features which can be used to guide motion correction in this direction. The relationship between image resolution and registration accuracy should be explored. Such information would be helpful for designing registration algorithms and interpreting registration accuracy within the context of the resolution of the dataset. Furthermore, computation time may be optimized by registering at a resolution that is known to permit a desired level of registration accuracy. Possible effects could be investigated using a set of unenhanced MR image volumes acquired at various resolutions in each coordinate direction. An identical sequence of motion could be simulated for each resolution level. A conventional registration algorithm would then be applied to each dynamic image series to correct for motion. Since the


60

motion characteristics are known, the accuracy of each set of registrations can be determined in each coordinate direction. The resolution of each dataset can then be compared to the associated registration accuracy in order to search for potential trade-offs.

3.2.3

Application and Testing

The proposed registration technique was developed and tested using abdominal DCEMRI data. Although the abdomen contains many different structures with unique motion and enhancement characteristics, the technique should be applied to other DCE-MRI sites (e.g. cardiac, breast, brain, prostate) to investigate broader applicability. The intensity correction presented in this work does not depend on the source of the contrast-enhanced data. Consequently it may be of use for registration of data from other modalities such as ultrasound or computed tomography. Work would need to be done to fine tune the algorithm for each modality particularly with respect to details such as regularizer weighting, number of optimizer iterations, and multi-resolution optimization specifications.

3.2.4

Effects of Registration and Motion on Analysis

Regardless of the source of the contrast-enhanced data, investigation should be performed into how registration affects contrast-agent dynamics. Recall that registration should remove the motion from the image that is caused by the patient while preserving the dynamics of the flowing contrast-agent. Since these dynamics are the focus of subsequent analyses, it is important to understand how registration and motion may affect them. Registration of a dataset with known simulated motion provides an indirect way of ensuring that registration accurately reproduces known ground-truth dynamics. If the registration is completely accurate then the contrast-agent dynamics in the registered dataset are also accurate. However, the sensitivity of the dynamics and subsequent modelling to inaccuracy is unknown. This sensitivity could be explored by simulating


61

varying levels of motion within a static enhancing dataset and investigating the effects on intensity-time profiles and pharmacokinetic modelling. Proper investigation of sensitivity presents a challenging task as it would be influenced by the type of motion simulated, the imaging site, the imaging modality, and the pharmacokinetic model. However, such an investigation would help to answer the question of how accurate registration needs to be within the contrast-enhanced setting.

3.3

Summary of Thesis Contributions

A novel approach for registration of contrast-enhanced images has been developed and tested using abdominal DCE-MRI data. The technique employs a motion-informed intensity correction which attempts to simulate contrast enhancement in the moving image such that it matches the level of enhancement in the reference. The intensity correction facilitates registration by re-establishing the assumption that differences between the two images are solely due to motion. Quantitative efficacy was demonstrated through registration of a DCE-MRI dataset containing a known amount of simulated motion. These results were put into context using a novel comparative approach that involved construction and registration of a control dataset. The control dataset contained anatomical and motion characteristics that were identical to the simulation but contained no enhancement. Since the ideal contrast-enhanced registration algorithm is unaffected by enhancement, registration of the control using conventional techniques provides insight into how such an algorithm would perform. Registration of unenhancing data is a much simpler task and therefore the associated registration accuracy sets a demanding standard for comparative purposes. The proposed registration approach was shown to perform well with respect to this standard. In addition, the proposed method was shown to be more accurate than a competing registration technique.


62

Potential clinical utility was shown through improvements to pharmacokinetic model fits to clinical DCE-MRI data. Motion artifacts in parameter maps were significantly reduced indicating the utility of registration within the context of voxel-based analyses. Improved fits obtained through registration support the recommendations of [36] and [48] for performing registration prior to DCE-MRI analysis. The technique is robust in that it does not rely upon manual ROI delineations, pharmacokinetic models, or assumptions about the motion characteristics of image data. The approach is also independent of the temporal sampling of the image data and is therefore suitable for pair wise registrations. Furthermore, it does not make use of any properties of the imaging modality, and therefore may be applicable within other contrast-enhanced imaging contexts (e.g. contrast-enhanced computed tomography, contrast-enhanced ultrasound).

Appendices

63

Appendix A Pharmacokinetic Model Fitting The fitting of a pharmacokinetic model to DCE-MRI data involves multiple steps which are discussed in detail below for a measured time-varying AIF version of the standard Tofts model. STEP 1: The concentration of the contrast agent in the blood plasma, Cp (t), is determined. An ROI within an artery that feeds the tissue of interest is selected. The average intensity of the pixels (voxels) in the ROI is calculated at each point in time during the scan. These values constitute the SPGR signal for the blood plasma, Sp (t). Knowing Sp (t), the flip angle α, the repetition time TR, and the equilibrium longitudinal relaxation magnetization Mp0 , the SPGR equation can be solved for the T1 of the blood plasma, T 1p (t). T 1p (t) =

Sp (t) − sin(α) Mp0 )]−1 −T R[log( Sp (t) − sin(α) Mp0 cos(α)

(A.1)

Mp0 can be determined for the ROI by performing some preliminary MRI experiments. The ROI can be imaged before the injection of contrast agent at several different flip angles. The SPGR equation can then be solved for T 1p0 , which is the relaxation time of the plasma when the contrast agent is not present. Using T 1p0 , the SPGR signal from the ROI from a single time point (before contrast injection), α, and TR, Mp0 can be 64

Appendix A. Pharmacokinetic Model Fitting

65

determined using the SPGR equation. Cp is then given by

Cp (t) =

1/T 1p (t) − 1/T 1p0 r(1 − Hct)

(A.2)

where r is the relaxivity of the contrast-agent which characterizes its ability to change the T1 relaxation times of neighbouring spins. If the ROI used to determine Cp (t) is not directly within the tissue of interest, a small time offset can be applied to Cp (t) to account for the time it takes for the contrast agent to flow from the ROI to the tissue. STEP 2: Initial values of the pharmacokinetic model parameters are guessed or estimated. For the standard Tofts model these parameters would correspond to K trans and ve . These values serve as a starting point for subsequent model fitting and optimization procedures.

STEP 3: The pharmacokinetic model is used to predict the concentration of the contrast agent in the tissue of interest. In the case of the standard Tofts model with a delta function AIF, equation 1.2 can be solved [59] for Ct to give Ct (t) =

K trans Cp0 (1 − e−kep t ) kep

(A.3)

STEP 4: The SPGR signal that would result from the predicted concentration of contrast agent in the tissue is calculated. In order to do this, T 1t (t) must first be calculated using Ct (t) and the expression T 1t (t) = [

1 + ve rCt (t)]−1 T 1t0

(A.4)

T 1t0 corresponds to the T1 of the tissue before contrast-enhancement and can be determined in a similar way as T 1p0 . Knowing T 1t (t), the SPGR signal for the tissue can be

Appendix A. Pharmacokinetic Model Fitting

66

predicted (Stpred ) using the SPGR equation. Stpred (t)

=

Mt0 (1 − e 1−e

−( TT1R ) t

−( TT1R ) t

)sin(α)

(A.5)

cos(α)

STEP 5: The actual SPGR signal from the tissue of interest in the DCE-MRI data, St , is compared to Stpred (t). A common approach is to use the chi-squared or reduced chi-squared score as a metric: 2

χ =

∑ (St − Stpred )2 i

χ2red =

σ2 1 χ2 DOF

(A.6)

(A.7)

where DOF is the degrees of freedom and is equal to the number of observations less the number of fitting parameters less 1. St is sampled in different ways depending on whether mean-value ROI or voxel-based analysis is being performed. STEP 6: The value of the metric is optimized by updating the estimates of the model parameters in an iterative optimization loop. Minimization of the chi-squared scores increases the agreement between the observed and model predicted SPGR signal for the tissue. If the model accurately predicts the true SPGR signal, then the associated model parameters are also assumed to be accurate (e.g. K trans , ve ).

Appendix B Principal Components Analysis and Imaging Principal components analysis (PCA) is a method of expressing data in terms of a coordinate system for which the basis vectors correspond to the variance components in the data. The first coordinate is in the direction of the maximal variance of the data, the second coordinate is in the direction of the second variance component, and so on [42]. By ordering the basis vectors (and the related components of the data in each direction) by the variance, the primary trends in the data can be isolated. Data compression and random noise reduction can also be achieved through reconstructing the data using only the principal (high variance) components. PCA techniques have previously been used during registration of DCE-MRI data [42]. The approach involves reinterpreting the image data by considering the intensity of each pixel or voxel as a time-dependent function. Let B represent a 2-dimensional DCE-MRI dataset and Bij (t) correspond to the pixel intensity at spatial index (i, j) and time index t = 1 . . . tmax . Now define ¯ij . Bij∗ (t) = Bij (t) − B

67

(B.1)

Appendix B. Principal Components Analysis and Imaging

68

Bij∗ (t) for t = 1 . . . tmax comprises a set of tmax vectors each with tmax elements. The tmax by tmax covariance matrix is calculated for this set of vectors. The eigenvectors of the covariance matrix are computed, and then sorted in order of decreasing eigenvalue. Denote the set of ordered eigenvectors by En , where n = 1 . . . tmax . B ∗ can be expressed as a linear combination of these eigenvectors.

Bij∗

=

t∑ max

cnij En

(B.2)

n=1

where cnij = Bij∗ · En . Time-dependent trends in the pixel intensities bias the covariance matrix. These trends are reflected by the En whose associated eigenvalues are high. Since the En are sorted in order of decreasing eigenvalue, the time-dependent trends are described by the first few eigenvectors. In contrast, temporally uncorrelated changes such as random noise or random patient motion, appear evenly throughout the covariance matrix. Consequently these changes are described by the later eigenvectors. Temporally uncorrelated changes can be omitted from the reconstructed dataaset by truncating equation B.2 for some β < tmax such that ∑

β