In-Plane Motion Correction in Reconstruction of non-Cartesian 3D-functional MRI

(1)

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

In-Plane Motion Correction in Reconstruction of

non-Cartesian 3D-functional MRI

Examensarbete utfört i Medicinsk Teknik vid Tekniska högskolan vid Linköpings universitet

av

Anette Karlsson

LiTH-ISY-EX--11/4480--SE

Linköping 2011

Department of Electrical Engineering Linköpings tekniska högskola Linköpings universitet Linköpings universitet SE-581 83 Linköping, Sweden 581 83 Linköping

(2)

(3)

In-Plane Motion Correction in Reconstruction of

non-Cartesian 3D-functional MRI

Examensarbete utfört i Medicinsk Teknik

vid Tekniska högskolan i Linköping

av

Anette Karlsson

LiTH-ISY-EX--11/4480--SE

Handledare: Maria Magnusson

Datorseende, ISY och radiofysik, IMH, Linköpings universitet

Olof Dahlqvist Leinhard

radiofysik, IMH, Linköpings universitet

Peter Lundberg

radiofysik, IMH, Linköpings universitet

Examinator: Maria Magnusson

Datorseende, ISY och radiofysik, IMH, Linköpings universitet

(4)

(5)

Avdelning, Institution Division, Department

Division of Computer Vision Department of Electrical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

Datum Date 2011-06-08 Språk Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport

URL för elektronisk version

http://www.control.isy.liu.se http://www.ep.liu.se ISBN — ISRN LiTH-ISY-EX--11/4480--SE Serietitel och serienummer Title of series, numbering

ISSN —

Titel Title

Korrigering av 2D-rörelser vid rekonstruktion av icke-kartesisk 3D funktionell MRI In-Plane Motion Correction in Reconstruction of non-Cartesian 3D-functional MRI

Författare Author

Anette Karlsson

Sammanfattning Abstract

When patients move during an MRI examination, severe artifacts arise in the reconstructed image and motion correction is therefore often desired. An in-plane motion correction algorithm suitable for PRESTO-CAN, a new 3D functional MRI method where sampling of k-space is radial in kx-direction and kz-direction and

Cartesian in ky-direction, was implemented in this thesis work.

Rotation and translation movements can be estimated and corrected for sepa-rately since the magnitude of the data is only affected by the rotation. The data were sampled in a radial pattern and the rotation was estimated by finding the translation in angular direction using circular correlation. Correlation was also used when finding the translation in x-direction and z-direction.

The motion correction algorithm was evaluated on computer simulated data, the motion was detected and corrected for, and this resulted in images with greatly reduced artifacts due to patient movements.

Nyckelord

Keywords MRI, fMRI, PRESTO-CAN, Motion Correction, Radial Sampling Pattern, Fourier domain, Circular Correlation, Gridding Reconstruction

(6)

(7)

Abstract

When patients move during an MRI examination, severe artifacts arise in the reconstructed image and motion correction is therefore often desired. An in-plane motion correction algorithm suitable for PRESTO-CAN, a new 3D functional MRI method where sampling of k-space is radial in kx-direction and kz-direction and Cartesian in ky-direction, was implemented in this thesis work.

Rotation and translation movements can be estimated and corrected for sepa-rately since the magnitude of the data is only affected by the rotation. The data were sampled in a radial pattern and the rotation was estimated by finding the translation in angular direction using circular correlation. Correlation was also used when finding the translation in x-direction and z-direction.

The motion correction algorithm was evaluated on computer simulated data, the motion was detected and corrected for, and this resulted in images with greatly reduced artifacts due to patient movements.

Sammanfattning

När patienter rör sig under en MRI-undersökning uppstår artefakter i den re-konstruerande bilden och därför är det önskvärt med rörelsekorrigering. En 2D-rörelsekorrigeringsalgoritm som är anpassad för PRESTO-CAN har tagits fram. PRESTO-CAN är en ny fMRI-metod för 3D där samplingen av k-rummet är radiell i (kx, kz)-planet och kartesisk i ky-riktningen.

Rotations- och translationsrörelser kan estimeras separat då magnituden av signalen bara påverkas av rotationsrörelser. Eftersom data är samplat radiellt kan rotationen estimeras genom att hitta translationen i vinkelled med hjälp av cir-kulär korrelation. Korrelation används även för att hitta translationen i i x- och z-riktningen. Test på simulerat data visar att rörelsekorrigeringsalgoritmen både detekterar och korrigerar för rörelser vilket leder till bilder med mycket mindre rörelseartefakter.

(8)

(9)

Acknowledgments

First, I would like to thank Maria Magnusson for the many hours spent helping and guiding me through this work. I would also like to thank Olof Dahlqvist Leinhard for all the thoughts and ideas along the way. A thanks goes to Peter Lundberg for giving me the opportunity to be a part of this field of research.

I would like to thank Anders Tisell and Gustav Ahlman for all the inputs during the coffee breaks. Finally, I thank my boyfriend Henrik Svensson for the time spent finding errors in this report and also for the support when the thesis work felt tough.

(10)

(11)

Introduction

At the Center for Medical Image Science and Visualization (CMIV), Linköping University and at the department of Radiation Physics at the University Hospital of Linköping, a new method of sampling k-space (Fourier domain) for functional Magnetic Resonance Imaging (fMRI) examination has been developed [1], [2], [3], [4], [5], [6]. The method is called PRESTO-CAN and as the name implies, it utilizes PRESTO (Principles of Echo Shifting with a Train of Observations) [7]. Instead of conventional sampling, the sampling is radial in the (kx,kz)-plane and Carteisan in the ky-direction. Also, the time resolution has been improved using golden angle ratio when sampling in the angular direction and combining that with a so called hourglass filtering technique. When sampling in a radial pattern, central k-space becomes denser sampled than outer k-space. With this property some samples in central k-space can be thrown away giving a better time resolution.

The implementation of a motion correction algorithm that works with PRESTO-CAN is presented in this thesis. The algorithm uses the radial sampling pattern but is not dependent upon the specific PRESTO-technique or the specifics of fMRI. This means that the algorithm could be implemented in other fields within MRI than fMRI with PRESTO.

1.1 Purpose and Goal

Motions during an MRI/fMRI scan can be devastating regarding image quality, see Figure 1.1. Therefore, in order to apply PRESTO-CAN in clinical examinations, an algorithm for reducing artifacts that occur when a patient moves is required.

The main goal of this thesis is to, within the reconstruction procedure, detect head movements occurring during a PRESTO-CAN fMRI examination and correct for those motions with the purpose to reconstruct an image with reduced in-plane motion artifacts.

(14)

1.2 Problem Formulation

There are six different kinds of rigid body movements (translation in x-, y- and z-direction and rotation around the x-, y- and z-axis). As mentioned before, an example on how movements can affect the image quality is shown in Figure 1.1.

(a) Without motion artifacts (b) With motion artifacts

Figure 1.1: A slice of a brain in the (x,y)-plane, illustrating how movements can affect the image quality.

This thesis will focus on in-plane motion correction in the (x,z)-plane which includes translation in the x-direction, translation in the z-direction and rotation around the y-axis. To reduce motion artifacts in the whole 3D volume is called through-plane motion correction, but the work required to implement this is be-yond the framework of this thesis. Even if in-plane motion correction is the main goal, the algorithm should not inhibit an extension to a through-plane motion correction algorithm in the future.

The sampling is performed in a radial pattern and is later gridded to a Cartesian coordinate system [8], [9]. The gridding method contains several steps:

1. Pre-weighting

2. Interpolation with a Kaiser Bessel filter to a Cartesian grid 3. Post-weighting

4. Compensation for filter kernel

The current pre-weighting is designed for an even sampling pattern in the an-gular direction, but when correcting for the estimated rotation, the radial sampling pattern does not apply to that condition any more. Therefore, in order to correct for rotation, the pre-weighting must be modified to handle uneven sampling as well.

The hourglass filter used in PRESTO-CAN is designed for an even sampling pattern and an alternative filter that uses the actual angle distance between the profiles is implemented, see Section 3.2.2.

(15)

1.3 Approach 3

1.3 Approach

Let the object/patient be located in the spatial domain. Then the Fourier domain of the object/patient corresponds to k-space. Translation in the spatial domain only affects the phase in the Fourier domain, and rotation in the spatial domain becomes a rotation in the Fourier domain as well.

These properties can be utilized to decouple the problem into four different steps. The first step includes estimating rotation by looking at the magnitude of the signal, which is unaffected by translation movements. Step two corrects for the estimated rotation in step one. Because correction of the rotation will cause an uneven sampling pattern, as mentioned in Section 1.2, step two has to include adjustments of the reconstruction algorithm, the gridding method. After step two there are only translation movements left. Step three estimates the translation and correction of the estimated translation is performed in step four.

1.4 Related Research

Due to the fact that rigid body movements during MRI/fMRI often result in images with poor image quality, a lot of research has been focused on motion correction for MR images in order to minimize those artifacts [10], [11], [12]. A few different ways of estimating affine movements are presented in this section. More examples on how to estimate movements can be seen in [10].

Some motion estimating techniques are based on using external sensors on the patient or by modifying the scan to include navigator echoes [13], [14]. Using external sensors/devices takes time, and to prepare the patient prior the fMRI scan may also be uncomfortable and perhaps even interfere with the response of the paradigm [12]. Because the goal in this thesis (Section 1.1) is to estimate and correct for movements during the reconstruction procedure, no more effort on external sensors and navigator echoes has been made.

Some methods that only use information from the original scan in order to reconstruct an image with less motion artifacts, are presented in [10], [11], [12]. One example is the DART algorithm described by L.C. Maas et al., where rotation and translation are handled as two separate problems [12]. The DART algorithm estimates rotation by first re-gridding the data points into a polar coordinate system, and then calculates the translation in angular direction by correlating two images and finding the maximum value which represent the movement. Correction of rotation is performed and the data points are re-gridded back to a Cartesian coordinate system. The translation is then, as rotation, estimated by correlation [12]. PRESTO-CAN uses sampling in a polar coordinate system, making the DART algorithm easy to implement within PRESTO-CAN.

Another sampling technique similar to the DART algorithm for estimating movements is called PROPELLER MRI and is described by J.G. Pipe [11], [15]. The sampling pattern of PROPELLER MRI is denser in central k-space and there-fore enables the opportunity to create several low-resolution images and comparing them with an average low-resolution image in order to reconstruct an MR image with reduced motion artifacts.

(16)

1.5 Report Outline

• In Chapter 2, background theory to the methods used in later chapters is presented.

• The implementation of the reconstruction method is described in Chapter 3. • In Chapter 4, the results from the implementation tests are presented. • The results from Chapter 4 are discussed in Chapter 5.

• Chapter 6 presents the overall conclusions and gives suggestions for future work.

(17)

Chapter 2

Background Theory

This chapter starts with an introduction to Magnetic Resonance Imaging (MRI) and gives a description of the PRESTO-CAN method. A few specific properties of the Fourier Domain are also presented.

2.1 Principles of Magnetic Resonance Imaging

This section gives an introduction of the field of magnetic resonance but note that it is far from a complete description. The text is based on [16] and [17].

2.1.1 Basic of Nuclear Magnetic Resonance

Protons, neutrons and electrons all have a physical property called spin. Atoms with unpaired nuclei give rise to a net nuclei spin greater than zero which also leads to a nonzero magnetic moment. The human brain consists mostly of fat tissue and water tissue, both constitutes the hydrogen atom (1H) which has a nonzero nuclei

spin. When a strong magnetic field is applied, the hydrogen atoms are aligned so that the magnetic moment is either parallel (lower energy level) or anti parallel (higher energy level) to the applied magnetic field. Slightly more atoms will be in the lower energy state giving a net magnetic moment parallel to the applied magnetic field.

The energy needed for the protons to transfer from the lower energy level to the higher energy level, ∆E, is described by,

∆E = hν, (2.1)

where h is the Planck’s constant (h = 6.626 · 10−34J s) and ν is the resonant

frequency of the spin, the so called Larmour frequency,

ν = B0· γ, (2.2)

(18)

where B0is the applied magnetic field and γ is the nuclear spin specific

gyromag-netic ratio. The gyromaggyromag-netic ratio for one hydrogen (a proton) is 42.58M Hz/T . When only the magnetic field B0 is applied, no signal will be detected since

all protons are in equilibrium. In order to detect a signal in the receiver coils, a magnetic field, B1 is applied perpendicular to B0. B1 is alternating with a

frequency similar to the Larmour frequency, making the protons to transition between the energy levels, giving a decrease of the net magnetization vector in the direction of the magnetic field, B0. The angle between the magnetization vector

and B0 is called the flip angle and it increases during the time B1is applied.

Repetition Time and Echo Time

When the B1field is turned off, an increased number of protons will transfer back

to the lower energy level, leading to a decreasing signal in the receiving coils. Before a new B1-pulse can be sent, the net magnetization must be at equilibrium

again. The time between the B1-pulses is called the repetition time (TR).

Right after the B1-pulse, the magnetization vector will be in phase. However,

due to the protons slightly different experienced magnetic fields, they will start to dephase, decreasing the strength of the received signal. After a certain time, t, after the first pulse, a new pulse can be applied, at 180◦, from the first pulse. This will, after the time 2t, also called the Echo Time (TE), make the net magnetization in phase again resulting in an echo signal detected by the coils.

Relaxation Times

In different tissues, the binding strength between the hydrogen atom and the other atoms differ. Therefore, the time for the protons to be back at the equilibrium net magnetization aligned with the B0field also differs. T1 is known as the time until

the net magnetization reaches 63 % of the equilibrium net magnetization. Also the time it takes to dephase will be different between different tissue types. In fact, there are two different types of decay in the transverse magnetization process. One is the pure proton interaction. The time it takes for the decay to dephase down to 37 % of the max transverse magnetization value is called T2. The other reason why the net magnetization in the (x,y)-plane will decay is due to inhomogeneity in the B0 field. The combination of these two decay processes is denoted T2*.

Spatial Encoding

If the magnetic field is constant all over the object, there is no way of knowing where in the brain the signal originates from. However, by applying gradient fields in the three spatial directions, this leads to spatially different magnetic fields and therefore also different spin frequencies. Knowing the history of the applied gradients, the spatially location of the signal received can be found. The signal that the receiver coil detects can be described by (2.3) and is equal to the Fourier transform of the spatial data which is also known as k-space in the field of MRI.

(19)

2.1 Principles of Magnetic Resonance Imaging 7

S(kx, ky, kz) = Z Z Z

s(x, y, z)e−j2π[kx,ky,kz]T[x,y,z]_dxdydz _(2.3)

The sampling of k-space can be done in many different ways by changing the way of how the gradient fields are applied and altered during the MRI examination.

2.1.2 Functional MRI

Functional MRI is a specific application within the field of MRI that maps which regions of the brain that are active when performing different tasks. One of the most commonly used methods to determine the active regions of the brain is the blood oxygen level dependent (BOLD) [16].

Blood Oxygen Level Dependent (BOLD)

When performing a task, more oxygen is needed in the part of the brain that is active. More oxygen is transported, even more than needed, so even if more oxygen is consumed, there will be a total increment of oxygen bounded to the hemoglobin existing in the active area. The ratio between oxygenated blood and deoxygenated blood increases and therefore changes the property of the magnetic moment of the blood.

There is a delay of the received signal because it takes a while for the human body to increase the blood flow after the patient starts performing the task. In order to detect changes in magnetic moment of the blood, long echo times are needed. In addition to that, the total scanning time needs to be relatively fast compared to MRI in order to detect the changes of activity over time. Therefore, in order to achieve better time resolution, fMRI images have poorer image quality, regarding image resolution, compared to conventional MRI images.

2.1.3 Example on Movements in functional MRI

In this section, a description of common movements in an fMRI-examination is given.

Translation Movements

When lying in an MRI-scanner, the patient is fixated in the y-direction, and in the x-direction the gravity helps keeping the patient in place. However, translation in the z-direction occurs more frequently and often in combination when performing the fMRI task [18]. Six different translation movement parameters of different patients can be seen in Figure 2.1.

Rotation Movements

The most common rotation movement is rotation around the y-axis, which often occurs when the patient concentrates on performing the fMRI task [18]. It is also

(20)

Figure 2.1: Six examples on how patients translate in the x-direction, y-direction and z-direction during an fMRI examination. Image source: [18].

quite usual that a patient tilts the head to the right or to the left, a movement seen as a slow increase or decrease lasting perhaps for several minutes. Some examples on rotation movements can be seen in Figure 2.2.

(21)

2.1 Principles of Magnetic Resonance Imaging 9

Figure 2.2: Six examples on how patients rotate around the x-axis, y-axis and z-axis during an fMRI examination. Image source: [18].

(22)

2.2 PRESTO-CAN

This section describes the theory behind a few specific properties within the PRESTO-CAN method, [1]-[6], that is important to know when implementing motion correction.

2.2.1 PRESTO

PRinciples of Echo-Shifting with a Train of Observations (PRESTO) is a Cartesian sampling technique suitable for fMRI and the BOLD effect that requires long echo times (TE). PRESTO has the possibility of acquiring TE > TR due to its special echo shift technique and provides T 2∗ weighted images in a fast way [7].

2.2.2 Sampling Pattern

PRESTO-CAN keeps Cartesian sampling in the y-direction while samples radially in the (x,z)-plane because radial sampling in k-space can provide both high spatial and high temporal resolution at the same time [1]. Figure 2.3 shows two different sampling patterns. The left is a so called Echo Planar Imaging (EPI) sampling pattern in Cartesian coordinates. One (kx, ky)-plane is sampled with a so called single shot, meaning that the whole plane is sampled before exciting a new pulse. For PRESTO-CAN, the sampling pattern is showed to the right in Figure 2.3. The 3D sequence in PRESTO-CAN can also be regarded as several 2D-slices stacked on each other in the y-direction. One such slice is shown in Figure 2.4a. The different planes (lines in 2D) are called profiles and are sampled from -kr to kr at different angular positions.

Figure 2.3: Two ways of sampling 3D data. The left shows a Cartesian sam-pling procedure and the right shows the PRESTO-CAN hybrid radial-Cartesian sampling pattern. Image source: [1].

When sampling k-space in PRESTO-CAN, the profiles are sampled with an angular increment close to the golden angle,

(23)

2.2 PRESTO-CAN 11

∆γorig= 180 ·√ 2

5 + 1 ≈ 111

◦_. _(2.4)

If sampled at the exact golden angle, the profiles will be distributed almost evenly around k-space, but will at the same time never visit the exact same angular position again. In order to be able to return to the same angular position again, the increment angle between two profiles is chosen to

∆γ = 360 N · M ≈ 360 N · ( N 2 · 2 √ 5 + 1), (2.5) where N and M are explained in the following.

N profiles is needed to sample k-space in order to create one image. The profile

numbers are given by n = 0, 1, 2, . . . , N -1. The angle representing profile number

n is n · 360/N . M then represents the integer value of the profile number, n,

corresponding to the angle closest to the golden angle. The increment angle, ∆γ, is then chosen to the angle closest to the golden angle, ∆γorig. The profile number

N +1 will start from the first angular position again, provided that the number of

profiles are chosen to be a prime number [19]. The angles are still evenly spread out, but over time several images are built on the same sample positions (but at different time positions).

In Figure 2.4, the sampling pattern is shown for N = 7 profiles. The number of M that gives the angle closest to γorig is equal to 2. The smallest angle difference between the profiles can be calculated by

∆Φ = 360

2 · N. (2.6)

With seven profiles as in the example, ∆Φ ≈ 26◦.

Note that, when sampling k-space, the data are sorted by time as shown in Figure 2.4b. However, before transforming to a Cartesian grid, see Section 2.2.4, the data have to be reordered according to the angle of the profiles, as shown in Figure 2.4c.

2.2.3 Hourglass Filtering

Better time resolution in the reconstructed images has been achieved by imple-menting a so called hourglass filter [2], [1]. When sampling with polar coordinates, the center of k-space is more densely sampled than the outer parts of k-space. For example, every profile samples the sample point at r = 0, but only one sample point is needed to represent the data. The best one to keep is therefore the profile at the time point of interest. An example of the result after hourglass filtering with seven profiles is presented in Figure 2.5. In this example, the time point of interest is profile n = 3. Starting with the oldest and the newest sample, which in Figure 2.5 corresponds to n = 0 and n = 6, a sample point is removed if the

(24)

(a) Overview of the sampled profiles in the (kx,kz)-plane

(b) Profiles sorted by time (c) Profiles sorted by the angle

Figure 2.4: Different ways to visualize the profiles sampled in the (kx,kz)-plane.

sampling theorem, (2.13), is fulfilled even without that point. Looking at the re-sult when the profiles are sorted by time, as in Figure 2.5b, the remaining sample points look like an hourglass.

(a) The Hourglass filter (b) Profiles sorted by time when filtered by the hourglass filter

Figure 2.5: The hourglass technique visualized in two different ways. Image source: [1].

(25)

2.2 PRESTO-CAN 13

Without making any compensation for the removed sample points, artifacts can appear during the reconstruction process due to the changes in sampling den-sity [1]. The removed sample points therefore get new values from the neighboring sample points using interpolating in the angular direction, see Figure 2.6.

Figure 2.6: Interpolation in the angular direction. Image source: [1]. The number of profiles for every image in a data set is always the same, and the newest sampled profile in image 2 is sampled at the same angle as the oldest profile in image 1. The angular distance between the profiles is equal along k-space. These properties make it possible to create the hourglass in advance and then apply the filter to the profiles corresponding to a certain image while the profiles still are sorted by time.

2.2.4 Gridding Reconstruction

If the data are sampled in Cartesian coordinates, the reconstruction can simply be performed using 3D inverse Discrete Fourier Transform (DFT), which is not the case when sampling in a non-Cartesian way. The reconstruction then either needs to transform the data to the spatial domain by a different procedure or to first re-sample the data into Cartesian grid and then use 3D inverse DFT. The advantage of DFT is that it is a fast procedure and therefore re-sampling the data and using inverse DFT is the fastest of the different alternatives [9].

There are different ways of re-sampling the data [9]. The one used in PRESTO-CAN is called the gridding technique, which works by interpolating the sample points to a Cartesian grid by using a small kernel. The gridding technique in PRESTO-CAN can be divided into four different steps, individually described below.

Pre-Weighting

Before interpolating to Cartesian sample points, where the distance between the sample points is equal in both x-direction and z-direction, weighting of the sample points is performed to compensate for the non-even radial sampling pattern.

(26)

K-space is divided into different segments based on Dirichlet tessellation, where segments are created depending on the distance between samples. Figure 2.7 illus-trates the Dirichlet tessellation where the gray areas represent segments belonging to a certain sample point. The weight given to a certain sample point is the inverse of its corresponding segment area.

Figure 2.7: An illustration of the Dirichlet tessellation used to determine the sampling density in PRESTO-CAN. The area of each segment is used in the pre-weighting function.

Due to even angular distance, the area of the segments at different radial positions can be determined by

W (Pr,φ) = ( r ·2π_N if r > 0, 0.52· π N if r = 0, (2.7)

where W (Pr,φ) is the weight given to sample point Pr,φ and N is the number of profiles (in Figure 2.7, N = 4).

Interpolation with Kaiser Bessel Filter to a Cartesian Grid

If the weighted data in polar coordinates are described as Fw(kr, kφ), the mapping to a Cartesian grid, Fw(kx, kz) is performed by convolution of the weighted polar data with a gridding kernel,

Fw(kx, kz) = C(kx, kz) ∗ Fw(kr, kφ), (2.8)

where C(kx, ky) is the Kaiser Bessel filter, which determines how each polar data point will distribute to the Cartesian grid. The kernel used in PRESTO CAN is a

(27)

2.3 Properties of the Fourier Domain (k-space) 15

4-point Kaiser-Bessel filter. The Fourier transform or apodization function of the Kaiser Bessel filter has the advantage of having relatively low side lobes compared to the main lobe, which makes it a good choice for avoiding aliasing in the image. In order to further avoid aliasing due to the side lobes, the sampling distance in x-direction and z-direction is chosen to be

∆kx, ∆kz = 1

2 · ∆kr, (2.9)

i.e. an oversampling of a factor two.

Post-Weighting

When using an ideal pre-weighting function, this step is actually not needed, since the purpose of the post-weighting is the same as for the pre-weighting, which is to compensate for inhomogeneities in sampling density. However, if the pre-weighting function is not ideal, post-weighting will be useful. Sometimes, the actual sampling density is hard to estimate, and then post-weighting also is useful [9].

The post-weighting works by calculating a normalization grid ρ(kx, kz), which registers all contributions from polar data points to Cartesian data points and is performed as

Fw(kx, kz) =

Fs(kx, kz)

ρ(kx, kz)

. (2.10)

After this, the 3D inverse DFT can be performed, giving fw(kx, kz).

Compensation for Filter Kernel

The last step includes compensating for the filter kernel in the re-sampling step. If this compensation is not performed, uneven intensity along the image will occur due to the appearance of the Fourier transform of the Kaiser Bessel filter. The compensation can be calculated by (2.11) since convolution in Fourier domain corresponds to multiplication in spatial domain, see Table 2.1.

f (kx, kz) =

fw(kz, kz)

c(x, z) (2.11)

2.3 Properties of the Fourier Domain (k-space)

As mentioned in Section 2.1, the signal is generated in the Fourier domain. Prior knowledge regarding the properties of the Fourier Transform is required to un-derstand this thesis work. In this section, important properties of the Fourier transform when reconstructing a PRESTO-CAN fMRI examination with reduced motion artifacts are described.

(28)

Spatial Domain DFT Domain (k-space) Translation: f (x − a, y − b) e−j2πN(akx+bky)_{· F (kx}_{, k}_y)

Rotation: f (cos(θ) − sin(θ)

sin(θ) cos(θ)  x y ) F (cos(θ) − sin(θ) sin(θ) cos(θ)  kx ky )

Added phase: e−j2πN(ax+by)· f (x, y) F (k_x− a, ky− b)

Correlation: (f ∗ g)(x, y) F (kx, ky) · G∗(kx, ky)

Multiplication: f (x, y) · g∗(x, y) (F ∗ G)(kx, ky)

Table 2.1: A few 2D Fourier Transform theorems.

2.3.1 Movements in Fourier Domain

The patient moves in the spatial domain, but the signal is measured in the Fourier domain. A patient that translates in the MR camera corresponds to adding a phase to the signal in the Fourier domain, while a patient rotating in the spatial domain corresponds to a rotation in the Fourier domain as well. In Table 2.1, the 2D Fourier transform theorems used in the thesis are presented.

2.3.2 Sampling in Fourier Domain

Here, the theory regarding sampling in the Fourier domain is presented.

The Sampling Theorem

In order to avoid aliasing, the sampling theorem must be fulfilled. If a certain field of view (FOV) of the patient is desired, the sampling distance must fulfill

∆kx,y,z≤ 1

F OVx,y,z

(2.12) when sampling Cartesian.

In PRESTO CAN however, the sampling is performed in polar Coordinates in two directions. In the radial direction, (2.12) still applies. The number of profiles to choose instead needs to fulfill

Nθ≥

π

2Nr, (2.13)

where Nθ is the number of angular positions, which is equal to the number of profiles and Nr is the number of sample points with the sample distance _{F OV}1

x in

(29)

2.3 Properties of the Fourier Domain (k-space) 17

The Frequency Components

Different regions in k-space correspond to different image components. Low fre-quency components (near k-space center) correspond to slow variations in the images, which are more frequently sampled in PRESTO-CAN than outer k-space. Outer k-space, however, corresponds to edges and details in the image. Rigid movements can be detected without the details corresponding to high frequency components. Therefore, a low pass filter combined with the hourglass filter im-proves the temporal resolution.

(30)

(31)

Chapter 3

Method

This chapter starts with a description of the motion correction algorithm, followed by a description of the modifications needed in the gridding algorithm and the hourglass filter when used for motion correction. The implementation of a phantom used for testing and validating the motion correction algorithm is also presented. Finally, the method used for testing the system, is described.

3.1 The Motion Correction Algorithm

This section starts with a description of the flow of the motion correction algorithm followed by detailed description for estimating movements and correcting for them. The flow chart in Figure 3.1 illustrates that the algorithm starts by estimat-ing rotation movements, see Section 3.1.1 for a detailed description of the rotation estimation. Rotation movements can be estimated without interference from trans-lation movements if only the magnitude of the data is considered since transtrans-lation only affects the phase of the signal, see Table 2.1.

In the first lap, a low pass filter can be applied with the purpose of achieving better temporal resolution for the translation estimation. The low pass filter sim-ply removes data points if the distance from k-space center is larger than a certain value. The low pass filter implemented removes every sample point sampled more than seven samples away from k-space center in radial direction. Applied com-bined with the hourglass filter, this leads to that information from only 13 profiles will help estimating the translation between two images.

However, before the hourglass can be applied, correction for the estimated rotation is performed, see Section 3.1.3. The reason why the correction of the rotation must occur before applying the hourglass filter is because the interpolation of the removed data points will be inaccurate if the profiles have rotated, see Section 3.2.2.

The hourglass can then be applied, followed by gridding reconstruction to Cartesian coordinates in the spatial domain. The translation is estimated on the images in the spatial domain, see Section 3.1.2.

(32)

Figure 3.1: A flow chart explaining the implemented motion correction algorithm.

After estimating the translation, the information gathered from the rotation and the translation estimation is used on the input data. The second lap starts by correcting for the translation movements on the input data, see Section 3.1.4.

The correction of rotation is performed again due to two reasons. If the low pass filter is activated during the first lap, the rotation correction must be per-formed again including all sample points. The second reason is the advantages on performing the translation correction before the rotation, see Section 3.1.4.

3.1.1 Estimating Rotation

The rotation between two data sets can be estimated by finding the translation in the angular direction between them [12]. In PRESTO-CAN, this is an advanta-geous method because the profiles are sampled in polar coordinates.

Before finding the translation in the angular direction, some modifications are made on the input data. First, the profiles are sorted in angular order instead of

(33)

3.1 The Motion Correction Algorithm 21

time order, compare Figure 2.4b with Figure 2.4c. This is done for every individual image. The absolute value of the data is calculated so that rotation is estimated only on the magnitude of the data. After that, some of the sample points are excluded from the estimation. Near k-space center, the signal is too similar in order to find any angular translation and further out in k-space the signal can be rather noisy [12]. Therefore, there is a trade-off between which samples in the radial direction to include in order getting the best result. A good first approach is to include every sample point fulfilling 1₄rmax < r < 3₄rmax where r is the radial value of the sample point and rmaxis largest r-value sampled, but tests are recommended in order to find the best estimation.

After these modifications of the data have been made, the rotation is esti-mated as illustrated in Figure 3.2. Applying inverse Fourier transforms, followed by conjugate and multiplication and finally one Fourier transform, corresponds to circular correlation of the non-rotated data (the reference image) and the ro-tated image. It is only the angular translation that is interesting, which is why the transform only needs to be performed in the angular direction. The angular translation between the two images is estimated by first summing in the radial direction (only the radial samples that were selected) and thereafter calculating the number of pixels between the midpoint and the position of the max point.

Figure 3.2: Illustration of the rotation estimation.

Sub Pixel Resolution using Parabolic Fit

To achieve a sub pixel resolution, parabolic fit can be applied. The max point together with the two points closest to the max point are adapted to a second degree curve. The max point on the curve gives the sub pixel estimation of the rotation.

(34)

3.1.2 Estimating Translation

As mentioned earlier, the hourglass filter and the low pass filter can be applied before estimating translation, which may result in a better time resolution. Trans-lation estimation is performed on the images in the spatial domain.

The method of estimating translation is similar to the method of estimating rotation. It is illustrated in Figure 3.3 and uses correlation between the first image (the reference image) and the translated image. The correlation theorem in Table 2.1 is utilized. The two images are 2D Fourier transformed, and the conjugate is applied to one of them before multiplication of the two Fourier transforms is performed. When transforming back to the spatial domain, the translation is estimated as the vector from the center of the image to the correlation peak.

Figure 3.3: Illustration of the translation estimation. Note that all images are presented by the real values of the data.

Sub Pixel Resolution using Parabolic Fit or Zero Padding

Achieving sub pixel translation estimation can be performed by either parabolic fit, as in rotation estimation, or by zero padding.

Adding zeros in the Fourier domain leads to interpolation of more sample points in spatial domain, giving a better spatial resolution. Zero padding is a method that is more robust than parabolic fit when interpolating in the spatial domain

(35)

3.1 The Motion Correction Algorithm 23

but requires a lot of computational power, which is why zero padding more than by a factor four is not to be recommended.

3.1.3 Correcting for Rotation Movements

The estimated rotation for every image compared to the reference image (measured in number of profiles) is placed in a vector. However, the interesting part is estimating which profiles that have rotated compared to each other.

Time Point of Interest for Different Images

In Figure 3.4, the profiles are sorted by time. The reference image, image 1, uses the profiles n = [0, ..., 6] and the time point of interest is determined to be n = 3. Image 2 uses the profiles n = [1, ..., 7] and the time point of interest is profile

n = 4.

Figure 3.4: An illustration of the time points of interest for two different images. If a rotation is detected between the first image and the second image in the ex-ample in Figure 3.4, it is considered that the patient has rotated between sampling profile three and sampling profile four.

For the first half of the profiles in the reference image, the patient is assumed to be still as is assumed also for the last half of the profiles corresponding of the last image. This assumption is made because no image has those profiles as the time point of interest.

Sort the Profiles in Correct Angle Order

(36)

Since rotation in the spatial domain corresponds to rotation in the Fourier domain, the profile has been sampled at a different angle than it was supposed to. The resulting real angle for the profile is calculated by adding the supposed angle to the estimated angle. After that, the profiles are ordered by the real angle increment.

Follow the Movement of the Patient

When looking at each profile separately, there is an interesting approach that can be implemented. Instead of correcting each image according to the reference image, the images can also be corrected so that they rotate with the patient.

For example, if the patient rotates between profile three and profile four in Figure 3.4, the profiles, n = 4, 5, 6 is rotated in the first image as usual. In the second profile, however, the profiles rotating are n = 1, 2, 3 instead of n = 4, 5, 6, 7 which would have been the case if every image is corrected according to the reference image.

3.1.4 Correcting for Translation Movements

The translation is, as rotation, estimated by comparing the images with the ref-erence image. In order to determine which profiles that have been translated in each image, a similar procedure as for rotation is performed. The result from the translation estimation is a vector in the (x,z)-plane.

When correcting for the translation movements, it is preferable to keep the data ordered by time and therefore correct for the translation movements before the rotation correction. This is because then no reordering of the estimation vectors is needed. However, when the estimation of the translation is performed, see Figure 3.1, the algorithm already has corrected for the rotation movements (which is obligated in order to estimate translation). The rotation is unaffected by in which the order the correction is made, but that is not the case for translation ((T rans) · (Rot) 6= (Rot) · (T rans)). In order to correct for the translation, the estimation vector needs to be modified as

mx mz = cos(φ) sin(φ) − sin(φ) cos(φ)  m0 x m0_z , (3.1)

where [mx, mz]T corresponds to the estimated translation used later in the trans-lation correction, φ is the estimated rotation for that profile, and [m0x, m0z]T is the translation vector estimated by the algorithm for that profile.

After that, the translation is corrected by applying a phase to the data in the Fourier domain according to

FT(kr, kφ) = e(−j2π(mxkx+mzkz))· F (kr, kφ), (3.2) where

(37)

3.2 Other Alterations 25 ( kx= krcos(kφ), kz= krsin(kφ). (3.3)

3.2 Other Alterations

As mentioned in Sections 1.2 and 3.1.3, when correcting for the rotation move-ments, the angular increment between the profiles is probably no longer the same as it would be if no movement occurred. In this section, the alterations needed in the PRESTO-CAN algorithm are described.

3.2.1 Alternating the Pre-Weight Function

The current pre-weight function is designed to handle different sampling density only in the radial direction. When correcting for rotation movements, the profiles no longer have equal sampling distance in the angular direction.

The alteration consists of adding an angular dependency to the pre-weighting function giving a certain sample point,

w(Pr,φ) = ( r ·φn+1−φn−1 2 if r > 0, 0.52_· π N if r = 0, (3.4)

where w(Pr,φ) is the weight given to sample point Pr,φ placed on profile n, N is the number of profiles, φn+1and φn−1 are the angle value for the profiles next to profile n. Equation (3.4) can be compared to (2.7) and the difference between the old and the new pre-weight functions is also illustrated in Figure 3.5.

(a) The old weight function (b) The new weight function

Figure 3.5: The old weight function and the new weight function used when one profile has rotated.

(38)

3.2.2 Alternating the Hourglass Filter

When implementing the hourglass filter described in Section 2.2.3, in combination with the rotation correction, a few problems resulting in poor image quality can occur. This problem and the solution of the problem are described in this section.

Problems with the Old Hourglass Filter

The angular distance between the profiles will most certainly not be the same over k-space any more. Therefore, when using the old hourglass filter, sample points may be thrown away so that the sampling theorem is not fulfilled any more.

When the hourglass application occurs before the rotation correction proce-dure, interpolating a new value to the removed sample point position will be from the supposed neighbor sample points. When later correcting for the rotation, the value given at a certain sample point is no longer built up by the new neighbors, but the old.

The problem is illustrated in Figure 3.6, where the patient rotates between profile five and profile six, which leads to that profile six is sampled between profile three and profile five instead of its supposed position. However, the new sample points interpolated to profile six uses the information from the supposed neighbors instead of its real neighbors.

Figure 3.6: An illustration of the problem that occurs when interpolating the new sample points before correcting for rotation to its actual place. The interpolation of the sample points are not given by the real neighbors.

Implementation of the New Hourglass Filter

To solve the problem with the difference in angular distance and the interpolation of new sample points, a few alterations are made.

The hourglass filter is implemented first after the rotation correction. For each image, the hourglass filter uses the real angular difference between the profiles

(39)

3.3 Simulated Data 27

when removing sample points. Looking at one sample point at a certain radial position, if the angular distance between the neighbors of that sample point ful-filles the sampling theorem, the sample point is removed. The hourglass filtering always starts with the profiles furthest away from the time point of interest, see Figure 2.5b.

One possible example is if a profile rotates and therefore becomes sampled at the exact angular position as another profile. Then all the sample points of the profile furthest away from the time point of interest is removed since that profile is not giving any new information.

The profiles are always sorted in the real angular order and the actual neigh-boring sample points in angular direction are used when interpolating the new points.

3.3 Simulated Data

The motion correction algorithm was developed using simulated data as input signal. This section describes how the input data was simulated within this thesis.

3.3.1 Simulating a Phantom

The signal from an fMRI scan is, as mentioned before, sampled in k-space. When simulating a phantom, a function whose Fourier transform is known is required. This is the case with a circular disc whose Fourier transform is the jinc-function, which is described by

jinc(r) = J (1, 2π · r)

r , (3.5)

where J is the Bessel functions of the first degree and r is the radial coordinate. A circular disc with its corresponding jinc-function in the Fourier domain can be seen in Figure 3.7. An ellipse is a scaled circular disc and its Fourier transform will be a scaled jinc-function.

The k-space phantom used in this thesis is built up by several jinc-functions simulated so that the corresponding ellipses in the spatial domain result in Fig-ure 3.8. Table 3.1 shows the parameter settings for the different ellipses that are building up the phantom. A slice of a real brain is seen in Figure 3.9a.

3.3.2 Simulating Movements

To be able to estimate rotation and translation, the simulation has to include movements of the profiles. These movements are also created in the Fourier do-main, applied to every jinc-function.

The translation is applied by adding a phase with a certain translation in the x-direction and the y-direction, using the theorem in Table 2.1.

(40)

(a) A small ellipse (b) A jinc-function

Figure 3.7: A small ellipse (3.7a) and its corresponding Fourier transform (3.7b)

Figure 3.8: The phantom reconstructed to the spatial domain

When adding rotation movements, it must be implemented in such a way so that the resulting ellipses in an image are rotated around the midpoint on the total image instead of around the midpoint of each ellipse.

3.3.3 Adding Phase to the Image

Due to inhomogeneities in the gradient fields there will be phase shift in the spa-tial domain. Therefore, the brain image becomes complex, even though the real brain certainly is real, see Figure 3.9. A simple phase shift, like in Figure 3.10, corresponds to a translation of k-space.

In order to get a more realistic simulation, a possibility to add phase to the spatial phantom is implemented. If the patient moves, the phase is considered to move in a similar way in order to not interfere too much with the estimation. The absolute value, the real part, and the imaginary part of the phantom where a small phase shift is added can be seen in Figure 3.10.

(41)

3.3 Simulated Data 29

Ellipse dx dz scalex scalez scale rot intensity Big ellipse: -5 -15 20 30 1 110 1 Small ellipse 1: 19 -15 13 20 1 40 1 Small ellipse 2: 21 13 13 10 1 30 1 Small ellipse 3: 4 7 9 25 0.6 120 1 Circle: -25 9 6 5 1 0 0.5 Thin line 1: 30 -15 5 1 0.8 20 -0.5 Thin line 2: 20 -27 5 1 0.8 45 -0.5 Thin line 3: 0 -32 5 1 0.8 90 -0.5 Thin line 4: -20 -25 5 1 0.8 135 -0.5 Thin line 5: -29 -15 5 1 0.8 160 -0.5 Inner object: 13 -20 11 13 1 0 -1

Table 3.1: The parameters used for the ellipses of the phantom in the spatial domain. The number of pixels relocated from the center corresponds to dx and dz. The parameters scalex and scalez are the size in respective direction of the ellipses where scale is the relative size between the ellipses. The parameter rot is the rotation of the ellipse given in degrees and intensity is the intensity value of the ellipse.

(a) Magnitude (b) Phase

Figure 3.9: The magnitude and phase of one slice of an MR image in the (x,z)-plane.

Figure 3.10: A visualization of how a translation in k-space affects the output signal in the spatial domain.

(42)

(43)

Chapter 4

Results

In this chapter, results from the implementation of the motion correction algo-rithm are presented. Note that every image is presented by its magnitude. First, the results of the motion correction algorithm with the recommended parameter settings is presented (Section 4.2), followed by the results from the different parts of the algorithm (Section 4.3).

Figure 3.8 is an illustration of the phantom without any motion applied. There-fore, that image is used as a reference image when analyzing the results presented within this chapter.

4.1 Common Parameters and Movements Used in

the Tests

The motion correction algorithm consists of different parts with many optional parameters. The parameter settings that are the same for all tests are listed in Table 4.1. The number of profiles in one image is chosen to the smallest prime number fulfilling (2.13).

Parameter

Profiles in one image: 127 Profiles in total: 254 Images in total: 128 Sample points in one profile: 160 Table 4.1: Parameter settings used in the tests.

(44)

Two different types of movements are applied when testing the algorithm:

Step movement: A movement is applied between profile 127 and profile 128.

This means that the first image and the last image consist of profiles that have not moved compared to each other. In image 64 and 65, however, half of the profiles have moved and the other half have not moved, making the applied movement harder to estimate. Figure 4.1 shows an illustration of the step movement. Observe that the size of the movement can change and that the movement looks the same for both rotation and translation.

Pulse movement: Here the simulated movement is applied for a short time and

then removed again. When a pulse rotation is applied, the profiles n = [128, ..., 140] (corresponding to image 65-76) has moved four profiles away from their supposed positions, see Figure 4.2a. When a pulse translation is applied, a linear increase of one pixel/profile is occurring from profile 125 to profile 129 (corresponding to images 62-66). The translation remains at four pixels translation until profile 133 (corresponding to image 70) and is followed by a linear decrease of one pixel/profile back to zero movement again, see Figure 4.2b.

(45)

4.2 Results From the Complete Motion Correction Algorithm 33

(a) Pulse rotation

(b) Pulse translation

Figure 4.2: An illustration of the short movements applied during the tests.

4.2 Results From the Complete Motion

Correc-tion Algorithm

The parameter settings used for the complete motion correction results are: • Sample points fulfilling 1

4rmax < r < 3

4rmax are included for rotation

esti-mation.

• A low pass filter and the hourglass filter are applied before estimating trans-lation.

• The output results are given both with and without activating the hourglass to the output.

(46)

• Parabolic fit is used to obtain sub-pixel resolution.

Figure 4.3 shows the result after applying a step movement (Figure 4.1) and letting the motion algorithm estimate and correct for the movement. Both uncor-rected and coruncor-rected images are shown.

(a) No motion correction

(b) Motion correction - no hourglass on output

(c) Motion correction - hourglass on output

Figure 4.3: The result before and after motion correction (with and without acti-vating the hourglass filter) where the applied movement is a step movement with four pixel translation in x-direction, four pixel translation in y-direction and a four pixel translation in the angular direction.

(47)

The corrected images have less motion artifacts than the non-corrected images. Image 128 in Figure 4.3a shows that the phantom is both translated and rotated. The corresponding corrected image in Figure 4.3b is registered to the same position as image 1. A few, so called, streak artifacts arise in the corrected images. When activating the hourglass on the output, the result for image 64 in Figure 4.3c gives a little poorer result than when not activating the hourglass.

Figure 4.4 and Figure 4.5 show the results after applying a pulse movement (Figure 4.2). The first (Figure 4.4) has both a rotation and a translation movement applied while the other (Figure 4.5) only has a translation movement applied.

The actual movement is hard to define looking at the non-corrected images, but artifacts are present in the image. The correction of the pulse movement applied in Figure 4.4 does not improve the image quality much, but a small improvement is achieved. For example, looking at image 1 in Figure 4.4a and comparing that to image 1 in Figure 4.4b, it can be seen that the streak artifacts are removed when using the motion correction algorithm. However, looking at the images when the hourglass is activated (Figure 4.4c) serious artifacts are shown in image 64 but the result for image 128 instead becomes much better. The reason for the fairly bad result in figure 4.4 depends on that the pulse rotation is not detected, see also Section 4.3.1.

The results from a pulse translation in Figure 4.5 give very good results re-garding image quality. The artifacts arising from the movement are almost entirely removed after the motion correction algorithm.

(48)

(c) Motion correction - hourglass on output

Figure 4.4: The result before and after motion correction where the applied move-ment is a pulse translation movemove-ment in the x-direction together with pulse rota-tion movement.

(49)

(c) Motion correction - with hourglass on output

Figure 4.5: The result before and after motion correction where the applied move-ment is a pulse movemove-ment with four pixel translation in z-direction.

(50)

4.3 Partial Results

This section presents the results from testing specific parts of the motion correction algorithm.

4.3.1 Estimating Rotation

Every sample point fulfilling 1₄rmax < r < 3₄rmax is included in the tests shown in Figure 4.6, Figure 4.7, and Figure 4.8. The estimation has been tested while alternating which radial positions to include. The results of estimationg a skip rotation is shown in Figure 4.6. It is discovered that other constitutions of radial positions do not improve the estimations. The estimation results for the pulse movement and the results from the step movement with parabolic fit in Figure 4.7 and Figure 4.8 are both worse than the result in Figure 4.6.

The step movement is tested by moving the profiles 1, 4, 10, and 20 profiles ahead corresponding to an angle movement of 0.0247, 0.989, 0.2474, and 0.4947 radians respectively. The result from all different step movements can be seen in Table 4.2, where the results are equal to the applied movement for all move-ments except when the profiles has rotated just one profile. In Figure 4.6 the step movement is four pixels.

A step movement of four profiles is used when testing the parabolic fit func-tion. Since the parabolic fit does not affect other things than giving a sub-pixel resolution, only one test is performed.

Figure 4.6: The estimation of the rotation with an applied step movement of profiles rotated four profiles.

The step estimation in Figure 4.6 is equal to the rotation applied. For the pulse estimation in Figure 4.7, no rotation is found. Using parabolic fit, the result becomes worse when the applied movement is an exact number of profiles, see Figure 4.8. However, the error is always less than 0.5 profiles when using parabolic fit for a step movement.

(51)

4.3 Partial Results 39

Figure 4.7: The estimation of the rotation with an applied pulse rotation move-ment.

Rotation (in profiles)

Image no: 1-61 62 63 64 65 66 67 68-128 Applied 0 0 0 0 1 1 1 1 Estimated 0 0 0 1 1 1 1 1 Applied 0 0 0 0 4 4 4 4 Estimated 0 0 0 0 4 4 4 4 Applied 0 0 0 0 10 10 10 10 Estimated 0 0 0 0 10 10 10 10 Applied 0 0 0 0 20 20 20 20 Estimated 0 0 0 0 20 20 20 20

Table 4.2: Applied step rotation movement of different sizes and their correspond-ing estimated rotation with slidcorrespond-ing window. No parabolic fit is applied.

(52)

Figure 4.8: Estimation of the rotation with an applied step rotation movement. Parabolic fit is activated.

(53)

4.3.2 Estimating Translation

The implementation of the translation estimation is the same in both x-direction and z-direction. Therefore, the tests in only one direction are presented. Ini-tial tests have been made validating that the estimations are not affected by the translation in the other direction.

In order to test if the hourglass filter, by itself or in combination with a low pass filter, can achieve a better time resolution, comparisons with the sliding window (meaning that no filter is applied) for both types of movements are performed. Figure 4.9 and Figure 4.10 show the estimation with the different parameters in a step movement and in a pulse movement respectively.

(a) Sliding window and Hourglass filter

(b) Hourglass filter + low pass filter and Hourglass filter

Figure 4.9: Translation estimation with a step movement of four pixels in the z-direction.

(54)

(a) Sliding window and Hourglass filter

(b) Hourglass filter + low pass and Hourglass filter

Figure 4.10: Translation estimation for a pulse translation in the z-direction.

combination with a low pass filter gives similar results while the sliding window results in a worse estimation. When the pulse movement is applied, the hourglass filter in combination with a low pass filter is superior to both sliding window and applying only the hourglass filter.

The estimation after applying the hourglass filter and the low pass filter is tested with a few different sizes of the step movement. The movements chosen are: 1, 2, 4, and 8 pixels. With larger movements, the phantom moves out of the FOV. The estimation result when applying different step translations are shown in Table 4.3. The estimation seems to be better for smaller translations than for larger.

A comparison between the parabolic fit and zero-padding is also performed. Zero-padding, as for parabolic fit, does not affect the estimation more than per-forming a sub pixel resolution. The estimation of activating parabolic fit or

(55)

acti-4.3 Partial Results 43

Translation (in pixels)

Image no: 1-61 62 63 64 65 66 67 68-128 Applied 0 0 0 0 1 1 1 1 Estimated 0 0 0 1 1 1 1 1 Applied 0 0 0 0 2 2 2 2 Estimated 0 0 1 1 1 2 2 2 Applied 0 0 0 0 4 4 4 4 Estimated 0 1 1 2 3 3 4 4 Applied 0 0 0 0 8 8 8 8 Estimated 0 1 1 3 5 6 7 8

Table 4.3: Applied step translation in the z-direction and the estimated transla-tion. Hourglass and low pass filtel activated. No sub-pixel estimatransla-tion.

vation zero-padding is shown in Figure 4.11. It can be seen that the results for parabolic fit and zero-padding are similar.

(56)

(a) Step movement of 3.75 pixels translation in the z-direction

(b) Pulse movement

Figure 4.11: Comparing parabolic fit with zero padding with two different kind of applied translation movements. Low pass filter and hourglass filter are activated.

(57)

4.3.3 Correcting for Rotation Movements

In Figure 4.12, the two different correction methods: standing still and follow the movement are shown. No estimation of the rotation is performed here; the applied movement is instead given directly to the rotation correction. No hourglass filter is activated here due to minimizing factors that can affect the image quality besides the rotation correction. Since the real applied movement is used for the estimation vector, the right correction will occur for every movement applied. Therefore, only one rotation, the step movement with a rotation of 10 profiles, is used when looking at the quality after rotation correction.

The rotation correction results in images with much better image quality. A few streak artifacts is however found in image 64, for both standing still reconstruction and follow the movement.

(58)

(a) No Motion Correction

(b) Motion Correction - Standing still

(c) Motion Correction - Follow the Movement

Figure 4.12: The result after applying the two different rotation correction al-ternatives for a step movement rotated 10 profiles in the angular direction. No estimations are included.

4.3.4 Correcting for Translation Movements

The result of translation correction of standing still and follow the movement can be seen in Figure 4.13. As for testing the rotation correction, no estimation is performed. The applied movement is directly given to the translation correction.

(59)

The applied movement is a step movement with a four pixel translation in the negative x-direction and a four pixel translation in positive z-direction. No artifacts at all are seen after translation correction.

(a) No Motion Correction

(b) Motion Correction - Standing still

(c) Motion Correction - Follow the Movement

Figure 4.13: The result after applying the two different translation correction alternatives for a step translation. No estimation is included.

(60)

4.3.5 The New Pre-Weighting Function

By applying a step rotation movement for a single ellipse, see Figure 3.7a, the old pre-weighting function is compared to the new pre-weighting function when looking at image 64. The Fourier transform of the ellipse, when no profiles have rotated is used as a reference image. If the pre-weighting is ideal, the result would look like the reference image. The effect of activating the post-weighting is also tested.

The result after applying the old weighting function and the new pre-weighting function is presented in Figure 4.14. The results after activating the post-weighting are also presented. The new pre-weighting function (Figure 4.14e) gives a smoother result than the old pre-weighting (Figure 4.14c), but still not as good as the reference image, Figure 4.14a. When also adding the post-weighting, no difference can be seen between the new pre-weighting, the old pre-weighting and the reference image (Figure 4.14f, 4.14d and 4.14b).

(61)

(a) No rotation, new pre-weighting, no post-weighting

(b) No rotation, new pre-weighting, post-weighting

(c) Rotation, old pre-weighting, no post-weighting

(d) Rotation, old pre-weighting, post-weighting

(e) Rotation, new pre-weighting, no post-weighting

(f) Rotation, new pre-weighting, post-weighting

Figure 4.14: The results of the image 64 (out of 128), for the k-space of a circular disc with a step rotation of five profiles in the angular direction using different weighting functions.

In-Plane Motion Correction in Reconstruction of non-Cartesian 3D-functional MRI

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

In-Plane Motion Correction in Reconstruction of

non-Cartesian 3D-functional MRI

In-Plane Motion Correction in Reconstruction of

non-Cartesian 3D-functional MRI

Examensarbete utfört i Medicinsk Teknik

vid Tekniska högskolan i Linköping

av

Abstract

Sammanfattning

Acknowledgments

Contents

Chapter 1

Introduction

1.1

Purpose and Goal

1.2

Problem Formulation

1.3

Approach

1.4

Related Research

1.5

Report Outline

Chapter 2

Background Theory

2.1

Principles of Magnetic Resonance Imaging

2.1.1

Basic of Nuclear Magnetic Resonance

2.1.2

Functional MRI

2.1.3

Example on Movements in functional MRI

2.2

PRESTO-CAN

2.2.1

PRESTO

2.2.2

Sampling Pattern

2.2.3

Hourglass Filtering

2.2.4

Gridding Reconstruction

2.3

Properties of the Fourier Domain (k-space)

2.3.1

Movements in Fourier Domain

2.3.2

Sampling in Fourier Domain

Chapter 3

Method

3.1

The Motion Correction Algorithm

3.1.1

Estimating Rotation

3.1.2

Estimating Translation

3.1.3

Correcting for Rotation Movements

3.1.4

Correcting for Translation Movements

3.2

Other Alterations

3.2.1

Alternating the Pre-Weight Function

3.2.2

Alternating the Hourglass Filter

3.3

Simulated Data

3.3.1

Simulating a Phantom

3.3.2

Simulating Movements

3.3.3

Adding Phase to the Image

Chapter 4