Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Investigating phase and artifacts in MRI : phase unwrapping, motional deghosting and fast dynamic flow… Chavez, Sofia Emilia 2000

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2000-565262.pdf [ 27.96MB ]
Metadata
JSON: 831-1.0085490.json
JSON-LD: 831-1.0085490-ld.json
RDF/XML (Pretty): 831-1.0085490-rdf.xml
RDF/JSON: 831-1.0085490-rdf.json
Turtle: 831-1.0085490-turtle.txt
N-Triples: 831-1.0085490-rdf-ntriples.txt
Original Record: 831-1.0085490-source.json
Full Text
831-1.0085490-fulltext.txt
Citation
831-1.0085490.ris

Full Text

I N V E S T I G A T I N G P H A S E A N D A R T I F A C T S I N M R I : P H A S E U N W R A P P I N G , M O T I O N A L D E G H O S T I N G A N D F A S T D Y N A M I C F L O W I M A G I N G By Sofia Emilia Chavez B. Sc. (Physics & Physiology) McGi l l University, 1991 M . Sc. (Physics) University of British-Columbia, 1995 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y in T H E FACULTY OF G R A D U A T E STUDIES D E P A R T M E N T OF PHYSICS A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA October 2000 © Sofia Emilia Chavez, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of ThySl'CS ^ / ^ S T W ^ O C W ^ V The University of British Columbia Vancouver, Canada Date Ochi 13 } 2-OGO DE-6 (2/88) Abstract In this thesis, the phase information contained in Magnetic Resonance (MR) images is studied while investigating the phase wrap and motional ghosting artifacts. The product of the work consists of novel processing techniques and data acquisition schemes. The work can be described in terms of three separate projects. The first project consists of studying the phase wrap artifact in MRI . A better character-isation of wrapped phase maps and the associated pole structures has led to the development of a new type of quality map, called a score map. Based on this, a new pole-guided-cutline method for phase unwrapping is realised. The second project consists of studying the motional ghosting artifact. An existing method of ghost suppression ( l pGEM) , relying on a gradient energy minimization (GEM) criterion, is improved. An analytical method of obtaining G E M is derived, leading to a new, more general method of ghost suppression: 2pGEM. The gradient energy of an image is proposed as a good quantitative evaluator of ghost suppression methods. The 2pGEM method of ghost suppression relies on the input of two images with signifi-cantly differing ghosts. A more efficient method, requiring less data acquisition, is presented. Although the proposed technique has not yet been implemented, due to some practical lim-itations, this work helps support the extended use of 2pGEM on ghosting artifacts caused by factors other than motion. The 2pGEM method is also tested on a new type of imaging sequence, called the OK-scan, which produces dissimilar ghosting patterns. The successful results emphasize the more versatile applicability of the 2pGEM method relative to the original, l p G E M . The last project, while in a preliminary stage, demonstrates how a dynamic flow profile can be extracted from the motional ghosts. Such profiles are clinically available by using cardiac gated or triggered procedures. The benefit of the proposed method is its speed and simplicity of implementation since no sophisticated hardware or patient preparation is necessary. Instead, the velocity-induced artifacts in ungated images are used to predict the sought profile. The results are promising since the profiles are comparable to those obtained using a clinical standard which relies on gating. ii Contents A b s t r a c t Table of Contents L is t of Figures L i s t of Tables Acknowledgments x 1 Introduct ion 1 1.1 Magnetic Resonance Imaging: a field map with magnitude and phase . . . . 1 1.1.1 A Magnetic Moment in a Magnetic Field 3 1.1.2 Relaxation Time Constants: T\ and 3 1.1.3 Detection of a Precessing Magnetization 5 1.1.4 Encoding Spatial Information in MRI 7 1.1.5 Imaging Echoes and Pulse Sequences 9 1.1.6 Description of /c-space 12 1.1.7 Image Characteristics: Contrasts, Resolution and Quality 14 1.1.8 Conclusion 15 1.2 Summary of Main Contributions and Thesis Outline 16 1.2.1 Summary of Main Contributions 16 1.2.2 Thesis Outline 18 2 A Study of the Phase W r a p A r t i f a c t i n M R I : A N e w M e t h o d of Phase U n w r a p p i n g 21 2.1 Understanding Phase Maps in MRI 21 2.1.1 The Nature of Phase Information 22 2.1.2 Phase Loops : Poles 22 2.1.3 Sources of Poles in MRI 27 iii 2.2 The Phase Unwrapping Problem 32 2.2.1 Review of the Phase Unwrapping Problem 32 2.2.2 Review of Phase Unwrapping Methods . . 33 2.3 The New Outline Phase Unwrapping Method 35 2.3.1 Method 35 2.3.2 Results 38 2.3.3 Study of Noise Effects 39 2.4 Three-Dimensional Phase Unwrapping 40 2.4.1 Limitations of Two-Dimensional Phase Unwrapping 40 2.4.2 Three-Dimensional Poles 43 2.5 Summary 44 3 The Use of Gradient Energy i n M o t i o n a l Ghos t ing Ar t i fac t Reduc t ion 45 3.1 Effects of Motion in MRI 45 3.2 Ghost Reduction Methods: Deghosting 48 3.3 Two-Point Deghosting for Superimposed Ghost Patterns 49 3.3.1 Introduction 49 3.3.2 Review of Existing Two-Point Interference Method 51 3.3.3 Analytical One-Parameter Minimization Technique: lpGEM 52 3.3.4 Comparative Results 54 3.3.5 Discussion of Gradient Energy Curves 54 3.3.6 New Generalized Two-Parameter Minimization Technique: 2pGEM . 57 3.3.7 Discussion on Window Size: Problem of Ghost Overlapping 60 3.4 A New Application for 2pGEM: Dissimilar Ghost Patterns 69 3.4.1 Introduction 69 3.4.2 Results 71 3.4.3 Magnitude Version of 2pGEM: General lpGEM 71 3.4.4 New Use of Gradient Energy: Conditional 2pGEM 74 3.5 Gradient Energy as a Measure of Image Quality 76 3.5.1 Introduction 76 3.5.2 Method 78 3.5.3 Results & Discussion 79 3.5.4 Conclusion 84 3.6 Ghost Suppression with One and a Half Excitations or Less 86 3.6.1 Motivation 86 3.6.2 Theory 86 iv 3.6.3 Method 88 3.6.4 Results 90 3.6.5 Discussion 92 3.7 Summary 97 4 Fast D y n a m i c F low Imaging: M a k i n g Use of the M o t i o n a l Ghos t ing Ar t i fac ts and Phase Information 99 4.1 Introduction 99 4.2 Theory 100 4.2.1 Review of Velocity Encoded MRI Phase Maps 100 4.2.2 Review of Dynamic Image Reconstruction: DIR 104 4.3 Method 105 4.3.1 Proposed Method for Dynamic Flow Imaging 105 4.3.2 Determination of the Ghost Spacing 108 4.3.3 Experimental Set-up : I l l 4.4 Practical Considerations 112 4.5 Preliminary Results 116 4.6 Discussion & Conclusion 117 5 Conclus ion 120 Glossary of Abbrevia t ions 124 Bib l iography 125 v List of Figures 1.1 Schematic diagram of the three-dimensional time evolution of the net magnetization vector, M , and relevant time dependent plots of the different components 6 1.2 Pulse Sequence for a slice selective gradient echo experiment 11 1.3 A typical 2DFT MRI pair of magnitude images 13 2.1 Phase map simulation with closed fringeline and no poles 23 2.2 Phase map simulation with open fringeline and a dipole structure. . . . . . . . . . . . 24 2.3 Figure demonstrating that the boundary conditions in MRI result in a phase map having the topology of a torus 25 2.4 Phase maps derived from a single random complex field demonstrating the changes in dipole connections with various phase offsets 26 2.5 Zoom of a random phase map showing several poles and their fringeline connections. . . . 28 2.6 Probability plot for the value of the sum of three wrapped phase differences, w y . . 29 2.7 Four phase maps, sampled with different resolutions and corresponding central horizontal profiles, demonstrate the creation of poles due to undersampling 30 2.8 Velocity encoded phase map of blood circulating in the heart with an open-ended fringeline and corresponding poles 31 2.9 Demonstration of streaks resulting from a conventional phase unwrapping algorithm. . . . 33 2.10 Four velocity encoded images and results from the pole-guided-cutline method 38 2.11 Velocity encoded phase map with unwrapped result, score map and relevant cuts 39 2.12 Simulations of a wrapped phase map with no cutlines and deteriorating signal-to-noise ratio. 40 2.13 Simulations of a wrapped phase map with a cutline along the bottom and deteriorating signal-to-noise ratio 41 2.14 Sketch representing the phase map of a fully wrapped vessel and corresponding three dimensional poles 42 2.15 Resulting pole rings in three dimensions showing the dipole connections 43 3.1 Comparison of magnitude images resulting from different implementations of the l p G E M . 54 vi 3.2 Close-up comparison of magnitude images resulting from different implementations of the l p G E M 55 3.3 Eg vs 9 curves in a region of noise 56 3.4 Eg vs 6 curves in a three regions of pure ghosts corresponding to three different ghost orders 57 3.5 Eg vs 8 curves in a region of 70 58 3.6 Vector representation of g\ and #2 with equal and unequal magnitudes 59 3.7 Comparison of l p G E M and 2pGEM on simulated images, 7i and 72, of superimposed ghost patterns with differing magnitude and phase 61 3.8 Comparison of l p G E M and 2pGEM on actual images, 7i and 72, of superimposed ghost patterns with differing magnitude and phase 62 3.9 The effect of different window sizes used to suppress flat ghosts for synthetic images. . . 64 3.10 Effect of varying window sizes on actual images 65 3.11 More images showing the effect of varying window sizes on actual images 66 3.12 Effect of varying window size on synthetic images showing ghost overlap 67 3.13 Effect of varying the window shape on synthetic images 68 3.14 Magnitude versions of two interleaved images with dissimilar ghosting patterns 70 3.15 Representation of two interleaved, fc-space trajectories for the OK-scan 72 3.16 Magnitude images showing the effectiveness of 2pGEM on images with dissimilar ghosting patterns 73 3.17 Magnitude results from applying the 2pGEM in complex form and in magnitude form to demonstrate the undesirable darkening of edges 75 3.18 Results of the Eg test and the conditional 2pGEM method 76 3.19 Images indicating the regions used for the regional Eg calculations 81 3.20 Magnitude images showing the results of applying several different processing techniques on the ghosted images of the first row of Fig.3.14 84 3.21 Magnitude images showing the results of applying several different processing techniques on the ghosted images of the second row of Fig.3.14 85 3.22 Figure demonstrating the undersampling ghosts that result from undersampling the k-space data and then Fourier transforming it back to image space 88 3.23 Figure representing three interleaved images, I\, I2 and 73, broken up into their various components: a common part, Jo, corresponding to the time-averaged signal, and three ghosts: gi, 52 and gz, respectively 89 3.24 Magnitude images of an intermediate image, 7 i e 3 o , and the corresponding fully sampled image, h 90 3.25 Figure demonstrating the canceling out of the undersampling ghosts of 7Q 91 vii 3.26 Comparison of the results obtained by applying the 2pGEM method on two fully sampled images (requiring 2 excitations) and applying it on a single fully sampled image with another intermediate image (requiring 1.5 excitations) 92 3.27 Figure demonstrating the problem of ghost overlapping for the efficient implementation of 2pGEM requiring only 1.5 excitations 93 3.28 Results demonstrating the implementation of the 2pGEM on the equivalent of only 1.25 excitations 96 4.1 Bipolar gradient waveform used to velocity encode MR images 101 4.2 Sketch of possible gradient waveforms for pulse sequences leading to velocity sensitized and velocity compensated phase images 103 4.3 Figure showing a magnitude ghosted image and the associated autocorrelation function with associated profiles 110 4.4 Figure showing two examples of a set of OK-sc&n interleaved flow sensitized and flow compensated images, with ghosts that interfere 115 4.5 Figure comparing dynamic flow profiles obtained by the proposed technique and using a standard clinical method 118 4.6 Figure comparing dynamic flow profiles obtained by the proposed technique in two cross-sections of the same hose 119 viii List of Tables 3.1 Table indicating various Eg values associated with previously shown images 80 3.2 Table of Eg values for quantitative evaluation of ghost suppression methods 83 ix Acknowledgments I would like to take this opportunity to thank the people that greatly influenced the work presented in this thesis. First and foremost, I would like to thank my supervisor, Dr. Qing-San Xiang for all his support. This thesis is the result of the many stimulating conversations and discussions I had with him throughout the past five years. He provided me with many ideas and was definitely a key player in my motivation to do the required work. I very much appreciate his enthusiasm and dedication to his work as well as the dedication, patience and respect he showed me. He has an incredible ability to simplify a complex problem and is also very knowledgeable. I feel privileged to have had the opportunity to spend many hours with him where all these qualities led me to gain insight into MRI and research in general, resulting in some of the most valuable learning experiences of my entire graduate career. I would also like to thank my husband Greg for his patience and encouragement since I know this thesis would not have been possible without his continual support. He was helpful both professionally and personally and for this I feel very lucky. His input greatly influenced this thesis and I know that the pride I feel in the outcome I owe to him. Finally, I would like to thank my son Mateo for his patience and for being such a good boy. He has done nothing but enrich my life with meaning, perspective, motivation and pride...all essential factors for the completion of this thesis. To my unborn baby, I have to thank you for being the final motivation I needed to finish, you are an inspiration! I would like to dedicate this thesis to my parents. I thank them for the education and guidance they offered me and the freedom they gave me to discover my own direction in life. x Chapter 1 Introduction Magnetic Resonance Imaging (MRI) originated in the 1970's when Lauterbur obtained the first image [1]. Competing against the well-established x-ray technology, well exploited since the turn of the 20th century, MRI is a relatively young imaging modality, clinically available only since the 1980's. The benefits of such a non-invasive and versatile modality have made it one of the most prominent diagnostic tools in hospitals today. The ultimate goal of the work presented in this thesis is to make better use of the information contained in an MRI image in order to improve its diagnostic ability. The hope is that this will lead to a better understanding of human physiology and disease, paving the path for improved prevention, treatment and cure of health-related issues. This introductory chapter is split into two parts. First, the basic principles of nuclear magnetic resonance (NMR) and MRI are given. This is meant to provide the reader, unfamil-iar with MRI, with sufficient knowledge of terminology, concepts and procedures pertinent to MRI. The references and brief discussions should be sufficient for an appreciation of the work contained in this thesis. For those already familiar with MRI, this section can be used as a reference of notations and definitions as well as a review. The second part of this chapter gives a summary of the main contributions contained in this thesis and emphasizes the importance and relevance of the work. It also outlines the structure of the thesis. 1.1 Magnetic Resonance Imaging: a field map with magnitude and phase MRI is based on the fundamentals of nuclear magnetic resonance (NMR) which has been around since the 1930's. NMR is a phenomenon occurring in magnetic systems where the nucleus of an atom possesses a magnetic moment due to its angular momentum. The magnetic moment can be viewed as a tiny bar magnet. The term resonance refers to 1 the fact that the magnetic moment will precess at a particular, resonant (or "natural") frequency, which is determined by the characteristics of the particular system. The first accurate measurement of a nuclear magnetic moment was performed in 1938 by I.I.Rabi and his colleagues [2, 3]. The existence of such a magnetic moment and the phenomenon of resonance is at the basis of all NMR and MRI theory and experiments. The simplest way to describe MRI is that it is an imaging modality which uses magnetic fields in order to affect some of the molecules in our bodies so that they will, in turn, respond by producing a signal that can then be detected. The most abundant molecules in our bodies that are imaged using MRI are those of water and fat. The primary focus of MRI is to image the unpaired protons of the hydrogen nuclei present in these molecules. The reason it is the hydrogen nuclei that are imaged is because NMR relies on the presence of a magnetic moment associated with an unpaired nucleon such as the proton in the hydrogen nucleus. The fact that there is a great abundance of water and fat in our bodies makes it so that these are capable of producing a signal greater than for other molecules in tissue. It is customary to refer to this as proton NMR which implies that it is the nuclei of the hydrogen atoms that are being imaged. There are three fundamental steps to the process of MR imaging. First, the protons that are to be imaged must be affected by an external transmitter system. This system is made up of both hardware, to provide the magnetic fields, and software, to control the application of the fields, in a time-varying manner that is carefully designed to cause expected and wanted changes in the state of the magnetic moments associated with the imaged protons. Such a process is what is referred to as encoding the MRI signal. The second step of MRI involves turning off the time-varying magnetic fields in order to allow for the now "excited" protons to produce an NMR signal. This signal must be received by appropriate hardware and constitutes the data acquisition. The third step of MRI is concerned with decoding the signal emitted by the protons and processing and displaying the images. This thesis is concerned primarily with this last step, involving signal decoding and image processing. However, in order to decode and process the data produced by an MRI scanner, it is essential to have a knowledge of the other two steps which created the data. In other words, the basic procedures of spatial encoding and data acquisition must be understood. In fact, most of the processing techniques presented in this work, depend on specific spatial encoding and data acquisition schemes which were developed for this purpose. In order to develop new techniques for signal decoding and image processing, the data must first be interpreted. For this, an understanding of the physics of NMR is required since it describes the effects of magnetic fields on the imaged protons as well as the signal emitted by them. The following is thus a brief review of some of the basics of NMR and 2 MRI. 1.1.1 A Magnetic Moment in a Magnetic Field A magnetic dipole moment, m, of a proton, is determined by its angular momentum, also called spin, S: m = 7S (1.1) where 7 is the gyromagnetic ratio which, for a proton, is equal to 26752 rad/(sec • gauss). (The boldfaced symbols represent vector quantities with regular typing representing their magnitudes.) In MRI, it is common to refer to the imaged protons as spins. When a magnetic dipole moment, m, is in the presence of a magnetic field, B , it has potential energy, Ep, given by: Ep = - m B . (1.2) Assume now that the static field B is in the z direction: B = Bzz. Quantum mechan-ically, the z component of S for a proton can only take on one of two values given by: Sz = ± h, where h = ^ and h is Planck's constant equal to 6.626 • 10~34 Joule • sec. There are hence two possible energy states for a proton in a static field Bz: Ep = ±7 Q ) hBz (1.3) where the lowest energy is associated with m being aligned with the field. The difference in the energy in the two states is therefore: AEP = jhBz. (1.4) Statistical mechanics says that m has a greater probability of being in the lower energy state (aligned with the field) as described by the so-called "Boltzmann's distribution". Therefore, exposing a sample containing a large number of spins to a static magnetic field, B , along the z direction (for enough time and with some energy absorbing medium) will eventually result in a net magnetization, M , along the z direction. This is the first step involved in obtaining an MR image and is called polarization of the spins. This net magnetization is a useful representation of the starting point of an MRI system, it can be used to explain a lot of what is observed once other time-varying magnetic fields are introduced. 1.1.2 Relaxation Time Constants: T\ and T 2 Before explaining how spatial information is encoded in MRI, the issue of relaxation should be addressed. There are two very important relaxation times or time constants that describe 3 the motion of spins exposed to magnetic fields. These time constants are characteristics of the spins being imaged and reflect their molecular environment. Relaxation processes happen to be at the core of the signal contrast present in MR images and an understanding is hence crucial for image interpretation. It was previously mentioned that given enough time and an energy-absorbing medium, a sample containing a large number of spins exposed to a static magnetic field, B , would become polarized, with a net magnetization, M , along the direction of B . The time-scale of the above process is characterised by the relaxation time, 7\. This time basically indicates the rate at which M , initially aligned with B along the z-axis, will decay back to this original position if tipped away from it. This process is driven by thermal equilibrium since M , aligned with B corresponds to the lowest possible energy state. Therefore, energy is required in order to tip M away from the z-axis and any process that causes 7\ relaxation must remove this energy from the system. The energy can be considered to be removed from the spins and put into the lattice (surroundings) hence this relaxation is referred to as the spin-lattice relaxation. Since the effect of this relaxation process is to realign M back along the z direction, also called the longitudinal direction, another term used for Ti is the longitudinal relaxation time. Other than the process of T\ relaxation, M will precess about B if it is tipped away from it. The rate of precession is given by the Larmor frequency, UIQ, which depends on the gyromagnetic ratio, 7, and the magnitude of the static magnetic field, B: coo = lB. (1.5) Once M is tipped away from the z-axis, it will have x and y components referred to as a transverse component: M T = (Mx,My). While Ti indicates the rate at which M reorients itself back to the z direction, T 2 indicates the rate at which M T shrinks or loses magnitude. Since M is made up of a sum of spins, its magnitude changes can be described in terms of the spins. Ignoring the effects of Ti relaxation, one could assume that tipping M away from the z direction causes a simple precession, where M T rotates about the z-axis with a constant magnitude. However, this will only be the case if all the spins making up M are coherently (in phase) precessing. Any factor that removes energy from the system and causes Ti, can also have an effect on the phase relationship between the spins. Any dephasing of these spins will result in a decrease in the amplitude of M T . Other mechanisms that do not remove energy from the system, such as interacting spins which mutually flip each other, can also produce dephasing of the spins hence T 2 is referred to as the spin-spin relaxation and it is usually shorter than Xi . Since the effect of the spin-spin relaxation process is the loss of the transverse component of M , another term for T2 is the transverse relaxation time. It is 4 important to make the distinction between dephasing processes due to spin interactions and those resulting from inhomogeneities in the main magnetic field. Both will result in a loss of M T and these are often lumped into a single effect characterised by the time constant T2*. However, the dephasing caused by main field inhomogeneities is recoverable while that caused by spin-spin interactions is not, as will be discussed in section 1.1.5. There is often ambiguity in this term and it is important to define which processes are taken into account when mentioning T<i-The phenomenological equation of motion describing the precession of M in a static magnetic field Bo and subject to the relaxation time constants, Ti and T 2, is called the Bloch equation. It can be written in terms of the three spatial components of M as: dMx Mx — = jMyBo-jr dMz Mz - M 0 dt T i where Mx, My and Mz are the components of M along the x, y and z axes respectively and M 0 is the magnitude of its longitudinal component at thermal equilibrium. The solution to these equations can be written as: ML(t) = M 0 + (M° - M0)e-^Tl MT(t) = M^e'^e^04 (1.7) where ML and MT indicate the longitudinal and transverse components of M , respectively, with their values at time t = 0 given by M° and My and where cu0 is the Larmor frequency. A schematic diagram of this time evolution of the net magnetization vector, M , and relevant time dependent plots are shown in Fig. 1.1. 1.1.3 Detection of a Precessing Magnetization While the above was concerned with describing how a magnetic dipole moment is affected by a magnetic field, the focus will now become how a magnetic dipole moment produces a signal that can be detected. The same principles of physics are invoked, only now there is a shift in starting point. Assume there is a net magnetization, M , precessing about the z-axis. This magnetization may have been tipped away from the z-axis by introducing a short-lasting strong magnetic field perpendicular to B 0 , in the x — y plane, causing a precession of M about its axis. Otherwise, as is the case in MRI, the tipping away from the z-axis can be obtained by applying a time-varying relatively weak magnetic field which 5 (a) (b) Figure 1.1: Schematic diagram of the three-dimensional time evolution of the net magnetization vector, M , in (b) and relevant time dependent plots of the different components in (a). Bo is the main magnetic field along the z direction, M 0 is the thermal equilibrium state of M and M° and M£ represent longitudinal and transverse components of M at time t = 0, respectively. (Figure taken from class notes of Q.S.Xiang [4] with permission.) is synchronized with the precessing M (i.e. at the Larmor frequency, CJQ, dependent on Bo)- For the main field strengths present in clinical MRI scanners (.5 Tesla -1.5 Tesla), this Larmor frequency is in the radio frequency range and hence this tipping, time-varying field is commonly referred to as the radio frequency field or RF pulse. Once M is tipped away from the z-axis by the RF pulse, it will have a transverse component, M T , varying in time as described by Eq.1.7. This means that if a conducting coil is placed with its axis aligned along either the x or y direction, a current will be produced in the coil due to the time-varying magnetic field through it. This phenomenon is described by Faraday's Law. The creation of current in a receiver coil by a time-varying magnetic field is the induction method of NMR which was first observed by Bloch et al. at Stanford University in 1946 [5, 6]. While Purcell et al. at Harvard University independently achieved NMR in condensed matter by using an absorption method [7], it is the NMR induction signal of Bloch et al. which is the basis of the modern MRI signal. Since this signal decays due to loss of M T caused by T 2 relaxation processes, it is commonly referred to as the free-induction decay signal (FID). Also, since the RF pulse is at the source of this signal, it is often referred to as the excitation pulse. Similarly, an excited magnetization will refer to M having a transverse component and hence producing a detectable signal. 6 It is important to note that the MRI signal is not only characterised by its amplitude or magnitude, but since it is precessing, it has a frequency associated with it. At any given instant in time, the MRI signal can thus be defined in terms of its magnitude and phase which indicate how strong it is and where it is in its cycle of precession, respectively. It can therefore be represented as a two-dimensional vector, with magnitude and direction, or as a complex number with magnitude and phase. MR images therefore spatially map this vectorial quantity and can be referred to as vector field maps. The phase is what distinguishes an MR image from that produced by other medical imaging modalities (such a x-rays and CT scans). It provides MRI with a unique capability for encoding a vast variety of clinically relevant information. 1.1.4 Encoding Spatial Information in MRI To make the distinction between the terms MRI and NMR is to note that MRI is an imaging procedure which means that the spatial distribution of some observable is mapped. More specifically, MRI uses the NMR signal produced by the protons in our bodies in order to determine their exact location. Other factors related to the chemistry of our bodies can also be deduced from their NMR signal. For the purpose of the following discussion, it will be assumed that our bodies are made up of stationary protons arranged in some particular distribution and the goal of MRI is to determine that distribution. In order to obtain an NMR signal, the sample or object being imaged must first be polarized as described previously. This is accomplished, in an MR imager, by a main magnet producing a strong, static magnetic field along the z direction of magnitude B0. This main field is meant to be homogeneous within the imaging region so that an entire object, placed in this region, experiences the same main field strength. The way that spatial encoding is introduced to the sample being imaged, is through the use of field gradients. These are magnetic fields which linearly change in strength along a particular direction. Although these gradients are in the x, y and z directions, this only indicates the direction of spatial variation. The magnetic field associated with each of these gradients is always along the same direction as the main field, the z direction. These field gradients, Gx, Gy and Gz, are superimposed on the main field and their effect is to add to or take away from the main field strength by a spatially dependent amount. This causes a predictable, spatial variation in the field experienced by the object being imaged. In practice, these field gradients are relatively very weak, having a maximum strength of 10 mT/m, creating a maximum change in the strength of the magnetic field that is typically three orders of magnitude less than the main field strength [8]. There is a variety of ways in which the three spatial directions can be encoded in MRI. 7 However, the most common and simple way to describe this is in terms of three distinct localization procedures. First, slice selection restricts the signal to a particular slice of the object. This is generally used to restrict the image to a single region along the z direction. This is accomplished by use of a field gradient along the z direction which then causes a spatial variation of the main field in that direction. Due to this field strength variation, the spins will have Larmor frequencies which also vary spatially. As already described, the RF pulse causes a tipping of the net magnetization into the transverse plane, referred to as exciting the magnetization. For the tipping to be effective, the RF must be matched to the Larmor frequency of the magnetization. If the Larmor frequencies vary spatially, the RF can be made to match a given range of Larmor frequencies (bandwidth) corresponding to a finite location along z, exciting only the magnetization at that location. This results in z localization commonly called slice selection. ; Now that Gz and the RF have allowed for the selection of a transverse plane, field gradients Gx and Gy are used to encode the spatial information within the plane. Two distinct procedures used for this are called phase encoding and frequency encoding. These can be explained in terms of Fourier imaging. Usually, x localization uses frequency encoding and y localization uses phase encoding although this is arbitrary and can be switched. The Fourier transform of a signal, made up of several frequencies, indicates the amount each frequency contributes to the signal and the relative phase relation between these. If an MRI signal is made up of spins with spatially varying Larmor frequencies, the amount of spins with each Larmor frequency can be found from the Fourier transform of the MRI signal. This is the basis of frequency encoding where a field gradient, Gx, is applied during the signal detection so that the spins in the slice precess at different frequencies, depending on their x position. A Fourier transform of the signal then indicates the amount of spins at each Larmor frequency which indicates the amount of. spins at each position along x. In other words, this gives a projection of all the spins present in the slice, onto the rr-axis. The phase information given by the Fourier transform is what allows for encoding in the other spatial direction. If a field gradient, Gy, is turned on for a short amount of time, ty, this will result in a spread, along y, of the Larmor frequencies of the spins for the short period, ty. The end result will be that the spins will have accumulated an amount of phase rotation dependent on their y position. By applying Gy for a period ty before the signal detection and then using Gx during signal detection for x localization, the projection of the spins within the slice, onto the x axis, will be affected. By repeating this procedure several times with varying Gy strengths, a collection of phase encoded x projections can be obtained. A two-dimensional Fourier transform (2DFT) of this two-dimensional collection of signals will result in x and y localization, producing a 2DFT MR image. 8 MR images are digital images, with the localization of the MR signals represented in a grid of small volume elements, called voxels for three-dimensional images, and pixels for two-dimensional images. Although the images typically represent a slice of tissue with a finite thickness (typically 1 to 10 mm), the term pixel is used most frequently. Since the MR image represents a vector field where each pixel can be assigned a complex number, any image acquired intrinsically has real and imaginary components which can be represented in two separate images. The intensity of each pixel in a given real or imaginary image is representative of the amplitude of the respective component of the MR signal emanating from that spatial position. The intensity is shown as a discretized shade of grey where the lighter the pixel, the more positive its associated value and the darker the pixel, the more negative its associated value. 1.1.5 Imaging Echoes and Pulse Sequences MR images are usually produced from a 2DFT of a series of signals. As previously described, for an NMR induction signal (FID), the MRI signal eventually decays. However, the MRI signal not only decays due to T2* relaxation processes, its decay is more complex since it depends on the transverse component of many spins precessing at a variety of frequencies due to the frequency encoding. As described for phase encoding, if a field gradient is turned on for a short period of time, the spins will accumulate an amount of phase rotation dependent on their position. This means that if they are initially in phase with one another, the application of a field gradient will cause a dephasing of the spins and this will have the same effect as the T2* processes: the magnitude of transverse net magnetization, and hence MRI signal, will decay. However, since this dephasing is controlled by the particular field gradient, a reversal of the polarity of the gradient will cause this process to reverse and produce a refocusing of the spins. Such a reappearance of the MRI signal is commonly referred to as the MRI echo and in particular, this refocusing by the field gradient is called a gradient echo. In practice, such a refocusing gradient pulse is commonly used after the slice selection Gz to undo the dephasing caused by the slice selection pulse. It is also used before the frequency encoding gradient, as a dephasor, so that the first portion of the frequency encoding gradient, applied during signal acquisition, serves to refocus the spins. The MRI gradient echo signal is thus usually collected over a period of time where the spins refocus and then dephase again such that it appears symmetric, with a central maximum magnitude corresponding to all the spins in phase. The frequency encoding direction is commonly referred to as the readout direction since the data acquisition or readout occurs during the application of the frequency encoding pulse. Another way to refocus the spins that are dephased due to field inhomogeneities is to 9 use the RF rather than the field gradient. This type of echo is called a spin echo. It was discovered in 1950 by Hahn [9] who observed that if he applied another RF pulse to a sample, a time r after the first excitation pulse, a second FID appeared a time r later. This discovery was one of the greatest contributions to magnetic resonance since the understanding of its formation led to a better understanding of the resonance phenomena in general and it was the origin of modern pulsed NMR experiments. Spin echoes are easier to understand and describe in terms of a 180° RF pulse used to trigger the echo. This is a variation of Hahn's original work and is the most common way that spin echoes are produced nowadays. For the following, only the dephasing due to main field inhomogeneities, causing a spatial dependent variation of the Larmor frequencies of the different spins, will be considered. The T 2 molecular level effects will be ignored since these are temporally random and the spin echo does not refocus any loss of signal due to that type of dephasing. The rotating axes are a useful representation in MRI and in particular, they facilitate the understanding of the spin echo phenomenon. Since the transverse magnetization is precessing about the z-axis due to the main field of strength B0, along that direction, x' and y' axes can be viewed as stationary in a frame rotating about the z direction at a rate given by u0 = jB0. The x' and y' axes are referred to as rotating frame axes. In practice, if the spins are excited by a 90° RF pulse applied along the rc'-axis, rotating the net magnetization, M , into the transverse plane along the y' direction, the individual spins start to dephase due to magnetic field inhomogeneities. If, after a time r, their phase is inverted by a 180° RF pulse about the y' axis and the dephasing process continues since their frequencies have not been altered, the spins will eventually refocus and be in phase again a time r after that. This sequence of RF pulses causes a spin echo at time 2r after the initial 90° RF excitation pulse and is denoted: 90°, — r — 180°/ — r — echo. This particular sequence of pulses is called the CPMG sequence named after Carr, Purcell, Meiboom and Gill. Other variations of this exist, also resulting in spin echoes but the CPMG sequence has many advantages [4] and is the most commonly used spin echo sequence. The spin echo has a greater amplitude than the gradient echo occurring at the same time since the spin echo refocuses many more processes (such as field inhomogeneities) than the gradient echo which only refocuses the dephasing caused by the gradient itself [8]. The sequence of steps involved in generating and collecting MRI data are referred to as pulse sequences since they are defined in terms of the application of magnetic field pulses: RF pulses and field gradient pulses. It is useful to represent the series of steps in a timing diagram which represents the timing of the different components of the sequence as well as their amplitude and duration. A typical pulse sequence is drawn as a series of lines, one for 10 each of the magnetic field pulses: RF, Gx, Gy and Gz. The RF pulse is usually symbolically drawn as a Gaussian or sine function with its direction and the associated tipping angle displayed next to it. The shape of the RF pulse is indicative of the slice profile. Since the steps represented in a pulse sequence are repeated several times for different phase encoding values, the phase encoding steps are usually all drawn simultaneously in such diagrams but it is understood that only one of the values for the pulse amplitude is used each time. The MRI signal and data acquisition window are also often indicated in such diagrams. The signal is usually drawn as an oscillating function that grows and decays due to the refocusing and dephasing of the spins as previously described. There are two characteristic times associated with every pulse sequence. One indicates the time between the initial 90° RF excitation pulse and the centre of the signal echo and is hence referred to as TE (echo time). The other characteristic time indicates the amount of time between consecutive 90° RF excitation pulses and is referred to as TR (repetition time). A typical pulse sequence for a gradient echo experiment is shown in Fig. 1.2. R F Signal slice selection phase encoding — T E -time frequency encoding (readout) Figure 1.2: Pulse Sequence for a slice selective gradient echo experiment. Each of the lines represents the timing, amplitude and duration of the magnetic field pulses: RF, Gx, Gy and Gz, involved in generating the MRI signal, as indicated. Time progresses from left to right. TE is represented by the distance between the centre of the RF pulse and the centre of the signal (echo). This sequence uses an RF excitation pulse applied along the rotating axis, x', and resulting in a tipping angle of 90°. 11 Each phase encoded induction signal is obtained by repeating the pulse sequence while changing only the amplitude of the phase encoding pulse. A simplification is to assume that the net magnetization is in thermal equilibrium (aligned along the main field) every time the pulse sequence begins. Since data collection occurs when the magnetization has been tipped in the transverse plane by the 90° RF excitation pulse, it is necessary to wait for Ti relaxation to occur before repeating the sequence and hence TR>> Ti. This causes the TR time to be much longer than TE, resulting in a lot of wait time. However, this ensures that every induction signal acquired corresponds to the exact same initial state and only the phase encoding is different. Other variations of this do exist in practice, where the tip angle is smaller than 90° or TR is shorter in order to achieve faster imaging. Variations of the TR affect the image contrast as will be discussed in section 1.1.7 The process of image interleaving makes use of the wait time occurring after data ac-quisition and before the next RF excitation pulse. Interleaving may involve using this time period to excite a different slice of the object in order to collect a phase encoded projection for a different image. It may also involve the application of a different pulse sequence on the same slice in order to obtain a different type of image of the same slice but providing different information. Even the same pulse sequence could be used to obtain a different image of the same slice due to the differing initial conditions. Furthermore, many images can be interleaved, depending on the exact TE and TR values of the sequences. Since most MRI involves multislice imaging, the procedure of interleaving is very commonly used. 1.1.6 Description of &-space Each of the data acquisitions, consisting of a specific phase encoding amount, have been referred to as phase encoded induction signals. Several such phase encoded signals are required in order to complete a full two-dimensional data acquisition which is required for x and y localization as previously discussed. This full two-dimensional data acquisition can be referred to as the data matrix. A 2DFT of this data matrix gives the two-dimensional MR image of the selected slice. Therefore, the data matrix is related to the image matrix via a 2DFT. A very useful representation of the data matrix relies on this relation and is called the &-space representation of the data matrix. The /c-space, like the MR image, is complex, with real and imaginary components which can be represented in two separate matrices or images. These are in fact the raw data matrices. The magnitude of these matrices is represented in Fig. 1.3 which shows a typical 2DFT pair. The dimensions of the matrices depend on the number of phase encode steps along the y direction and the number of sampling points taken during signal detection along the x direction. This gives the image resolution as discussed in the next section. 12 (a) (b) Figure 1.3: A typical 2DFT MRI pair of magnitude images, (a) shows the magnitude of the raw data matrices, corresponding to fc-space, where every row is a phase encoded induction signal from the entire slice, (b) shows the magnitude image corresponding to a 2DFT of the raw data matrix. (Figure taken from class notes of Q.S.Xiang [4] with permission.) The A;-space data, written S(kx, ky) is usually collected line by line, where each line is a phase encoded induction signal corresponding to a specific ky value. In this fc-space representation, the coordinates are related to the field gradients along x and y and time, Every time Gx is turned on for a period of time, the /c-space is traversed in the kx direction. When a phase encoding step, Gy is turned on for a brief moment, a shift along ky takes place in A;-space. The velocity of these moves in A;-space depends on the direction and magnitude of the gradients. Once the slice selection portion of the pulse sequence has been performed, the effects of the rest of the applied field gradients can be viewed as changing location in /c-space. During data acquisition, the path traced out by the readout gradient in A;-space is referred to as a trajectory. This representation is a very useful tool for interpreting the variety of pulse sequences available since they all can be translated into a set of trajectories, covering the /c-space. This representation can be used to describe the pulse sequence shown in Fig. 1.2. Drawing all the trajectories for the readout gradient corresponding to Fig.1.2 will result in a series of horizontal, parallel lines due to the following. Immediately following the slice selection by: kx = lGxt ky = jGyt. (1.8) 13 pulses, the system is at the central position in A;-space since no gradients in x or y have been turned on. The phase encoding gradient then corresponds to a translation in A;-space to a specific ky position, still at kx = 0. The frequency encoding readout gradient along x corresponds to a translation along kx at that ky position. The negative lobe preceding the gradient echo readout pulse causes a move to the far left of /c-space. The readout period then traverses &-space from left to right, tracing out a trajectory. When the next slice selection occurs, if one assumes the system has returned to the same initial conditions, one can return to the central position of ft-space and repeat the above steps, this time, starting with a move to a different position along ky. RF pulses used after the slice selection pulse, such as the refocusing 180° pulse of a CPMG sequence, also cause a move in &-space. For example, the spin refocusing 180° RF pulse along the y' direction corresponds to a inversion of the phase of the magnetization with respect to this axis. This translates in A>space to a jump from the position (kx,ky) to the position (—kx, —ky) since the A>space coordinate can be shown to be directly proportional to the phase of the magnetization. Al l these moves in A>space make the fc-space representation an essential tool for pulse sequence development and interpretation. 1.1.7 Image Characteristics: Contrasts, Resolution and Quality-It has already been noted that the same part of the body can produce MR images with different contrasts or information due to the application of a different pulse sequence. This is mainly due to the effects of T\ and T 2 relaxation processes on the intensity of the signal. Different contrasts can thus be obtained from the same pulse sequence by varying the characteristic times, TE and TR, of the sequence. In general, the longer the TR, the more longitudinal relaxation will be available for the next RF excitation pulse and hence the stronger the signal. In contrast, the shorter the TE, the less signal dephasing due to T 2 processes and hence the stronger the signal. If TR is chosen so that TR » Ti and TE is chosen so that TE << T 2 then all the protons in a given slice can be thought of as contributing fully to the signal and hence such an image represents the proton density. However, since Ti and T 2 vary for the many different tissues in the body (where the protons are in different molecular environments), the contrast between these tissues may be more or less emphasized by causing the signal to depend more or less on these relaxation processes by adjusting the values of TE and TR accordingly. T2-weighted images are produced by choosing TR » Ti, so that the contrast is independent of the Ti of the tissues and by choosing TE ~ T 2 of the tissues so that the contrast is highly dependent on T 2. On the other hand, Ti-weighted images are produced by choosing TE « T 2 and TR « Ti of the tissues so that the contrast is independent of the T 2 of the tissues but highly dependent on 14 their Ti. The spatial resolution of an MR image is determined by the field-of-view (FOV) and the number of samples taken in each direction. Since the spatial encoding is accomplished by the use of linear field gradients and discrete sampling in time, there is only a given distance over which the spins can be made distinguishable. This distance is called the FOV and in each direction, it depends on the corresponding field gradients, Gx and Gy, and the sampling time steps, tx and ty: FOVx 2?r T ' Gx • tx F 0 V - = T ^ V <1J» Any spin occurring outside the FOV will be indistinguishable from another spin, within the FOV, and will hence be incorrectly localized due to a process called aliasing. For Ny phase encode steps and Nx frequency encoded readout sampling points, the spatial resolutions in x and y are given by: Ax = FOVx N, X Ay = ™ i ( 1 , 0 ) The signal-to-noise ratio (SNR) is an important measure of image quality in MRI. It expresses the strength of the relevant NMR signal relative to the noise contaminating the image. It is usually defined in terms of the mean intensity, within a region-of-interest (ROI) on the tissue, divided by the standard deviation of the intensity within an ROI containing only background noise. Many factors affect the SNR in MRI by either affecting the amount of contaminating noise or affecting the average intensity of the signal. Some of these include the main field strength, the coil size, TR, TE, the slice thickness, the number of signals averaged (also called NSA) and the FOV [4]. 1.1.8 Conclusion Many basic concepts, procedures and terminology used in MRI have now been defined. It is hoped that this will provide the unfamiliar reader with enough information so that the content of the rest of the thesis can be appreciated. Many issues and concepts have been greatly simplified and others ignored. This is not a reflection of their importance in MRI and it is not implied that these issues and ideas are any less interesting than the ones discussed. The choice of content for this chapter was solely based on what was judged essential for further reading of this thesis. For those already familiar with MRI, this section can be used as a reference of notations and definitions as well as a review. 15 1.2 Summary of Main Contributions and Thesis Out-line 1.2.1 Summary of Main Contributions This section describes the main outcomes of the work performed for this thesis. The purpose is to emphasize the importance and relevance of the work and place it within a broader context. The main contributions of this thesis consist of novel techniques for processing MRI data with artifacts. In particular, the artifacts studied are the phase wrap artifact and the motional ghosting artifact which both can stem from motion occurring during an MRI scan. These will be further defined and explained in the following chapters. The image processing techniques presented in this thesis are the outcome of studying the nature and origin of these artifacts. Therefore, although they can be described as simple tools, their development and use relies on a deeper understanding of the physics of MRI. It is hoped that, by presenting the motivation and derivation of the techniques in..this thesis, a better understanding of the underlying physical principles of MRI can be gained, paving the path for new discoveries and ideas. The work performed can be classified into three different projects with the common focus of MRI artifacts. An image artifact is the term used to describe any corruption of the ideal image. There are many sources of artifacts in MRI [10] and in general, a thorough study of the origin of the artifact can reveal many important things about the physics involved in the imaging procedure as well as the object being imaged. Two of the three projects have the goal of eliminating an MRI artifact while the third project involves making better use of the information contained in the artifact. Another important commonality between the three projects is that they exploit the complex nature of the MR image. Either explicitly, for the phase unwrapping studies, or implicitly, for the motional deghosting and dynamic flow studies, the work involved in all three projects is based on the fact that the MR image carries important information in the phase of the signal. This makes the work of this thesis quite specialized in that the processing techniques were developed specifically for MR images. However, any other medical imaging modality or field of study where similar artifacts corrupt the data, impeding proper extraction of relevant information, could make use of the ideas of this thesis. This will be pointed out in the following, every time an idea with potential use elsewhere is presented. The first project studies the phase wrap artifact in MRI. The phase of an MR image carries very important information but this artifact has limited its potential use in clinical settings. In particular, the behaviour and origin of the phase inconsistencies called phase 16 loops or "poles" appearing in MR images, was explored. This led to a new formulation of the problem of phase unwrapping which deals with overcoming this artifact. The final product was the realisation of a new technique for phase unwrapping which can hopefully be used to make better use of the phase of MR images. This new phase unwrapping technique relies on the creation of a new type of image quality map [11] called a score map. Both the score map and the phase unwrapping technique proposed here have the potential of being useful in other fields of study where the imaging of phase is important (such as for radar interferometry, see Chapter 2). This new phase unwrapping technique was first presented at the Fifth Scientific Meeting and Exhibition of the International Society of Magnetic Resonance in Medicine (ISMRM) in 1997 [12]. An associated publication has resulted which contains most of the contents of Chapter 2 [13]. The second project studies the motional ghosting artifact which is a very common ar-tifact in MR images. The corruption caused by this artifact may be significant enough to compromise the diagnostic ability of MRI. There are three main contributions that resulted from the work for this project. First, an analytical expression was found in order to solve a minimization problem that was the key to a previously proposed method of ghost sup-pression for MRI. This resulted in a faster and more exact implementation of an already existing ghost suppression method. This implementation was first presented at the Sixth Scientific Meeting and Exhibition of the ISMRM in 1998 [14]. The second contribution was a new formulation of the problem which led to a more general form of the ghost suppression technique, extending its applicability and motivating the creation of a new way of obtaining ghosted images with dissimilar ghosting patterns. This new ghost suppression technique was first presented at the Seventh Scientific Meeting and Exhibition of the ISMRM in 1999 [15]. The third contribution comes from a more thorough study of the scalar quantity called the gradient energy which was extensively exploited by the aforementioned ghost suppres-sion techniques. A systematic study of the value of this measure on many different images with varying levels of artifacts resulted in a better appreciation of the information it can convey. This led to an expansion of its potential use since several new applications which could benefit from such a measure were explored. This is an important contribution since the gradient energy of an image is defined for magnitude images as well as those of complex nature. Hence, the ideas developed in this study could potentially be applied to many other areas of imaging. The gradient energy as a good measure of image quality was first proposed at the Eighth Scientific Meeting and Exhibition of the ISMRM in 2000 [16]. A possible limitation of the motional ghost suppression methods studied and developed is that they rely on the acquisition of two full images, in order to obtain a single result. There was therefore some motivation to make the procedure more time efficient by reducing the 17 amount of necessary data. This was motivated by the fact that some of the data in the two images seemed redundant. It was thought that the gradient energy could provide enough complimentary information to recreate equivalent results with the acquisition of less data. A nice idea was developed and it is presented in this thesis. This work was first presented at the Seventh Scientific Meeting and Exhibition of the ISMRM in 1999 [17]. Unfortunately, some limiting factors are yet to be overcome, if this idea is to become practical. The third project studies the potential of using the information contained within the mo-tional ghosting artifact, along with the phase information contained in MR images, in order to extract a dynamic flow profile in a fast way. This profile describes the time evolution of the average velocity of the signal in a given region of an image. This study is in its prelimi-nary stage with only few results and many unanswered questions. The results are, however, very promising considering the poor level of sophistication of the current implementation of the technique, where many factors have yet to be optimized and completely understood. The main contribution from the results of this study is that it was demonstrated that a new, very simple and fast method of dynamic flow imaging is possible and that the idea merits further investigation. The importance of this study is that it emphasizes that an understanding of the source of artifacts present in MRI can lead to the extraction of useful information from the images. Furthermore, it ties together the other two seemingly unre-lated projects in that it not only explores the motional ghosting artifact but relies on a good appreciation of the phase information associated with an MR image. 1.2.2 Thesis Outline This section describes the layout of the contents of the thesis which can be separated into three distinct projects as described above. The first project involves the study of the phase wrap artifact in MRI. The second project addresses the motional ghosting artifact and the third project makes use of this artifact, in combination with velocity encoded phase maps, in order to recreate a dynamic flow profile. Each one of these projects makes up a chapter of the thesis, the contents of which are further described below. Chapter 2 consists of the work performed for the phase wrap artifact studies. It provides an introduction to the nature of two-dimensional phase maps in MRI: a description of some important features of the behaviour of phase maps is given. It then proceeds to develop a new formulation for the problem of phase unwrapping. This is followed by a review of the current status of phase unwrapping methods in other fields of study as well as MRI. A new cutline method for phase unwrapping MR images is then presented with some typical results demonstrating its capabilities. This chapter "wraps up" with a discussion on a new promising direction for the study of phase unwrapping MR images relying on three 18 dimensions. This is presented as a possible extension of the proposed two-dimensional method and some simulations are given to motivate further investigation in this area. A manuscript, which has been reviewed and re-submitted in final form for publication [13], contains most of the contents of this chapter . Chapter 3 consists of the work performed for the motional ghosting artifact studies. It begins, in 3.1, with an introduction to the effects of motion on the MR image. This is meant to expand the ideas presented in the first section of this chapter which describe the basics of MRI for the ideal situation of a stationary object being imaged. It is then followed, in 3.2, by a review of some of the present ghost suppression methods in order to put the work into context. The following sections of the chapter: 3.3-3.6, are dedicated to the new techniques being proposed. They are each structured more or less in the familiar format: theory, method, results and discussion. The first two contributions mentioned above consist of the development of an analytical implementation to an already existing method of ghost suppression as well as a more general method which allows for less restrictions on the two ghosted input images required. These two techniques are developed and used in sections 3.3 and 3.4 which have been structured according to the relation between the ghosting patterns of the input images: either superimposed in section 3.3 or dissimilar in section 3.4. Section 3.5 then presents the work involved in studying the gradient energy as a measure of image quality. Since the gradient energy has been extensively used for the ghost suppression techniques, this section contains a quantitative analysis of the results presented in 3.3 and 3.4. The last section, 3.6, presents the idea of using less data in order to obtain similar ghost suppression with improved time efficiency. Chapter 4 consists of the work performed for the dynamic flow imaging studies. It begins by continuing to discuss other effects that motion, in particular blood flow, will have on an MR image. Section 4.1 builds on the basics of MRI described in the first part of the present chapter and in the introduction to Chapter 3. Section 4.2 provides the theory behind velocity encoded MRI phase maps as well as an already existing dynamic image reconstruction technique. This is followed, in 4.3, by a description of the proposed method and in 4.4, by a discussion on some practical considerations that merit further investigation. Sections 4.5 and 4.6 contain the preliminary results and a discussion, respectively. The final chapter of this thesis contains a general conclusion based on the outcome of the three projects undertaken. As is usual in research, the investigation of a given problem revealed several other interesting problems. Not only have some useful new techniques for processing MR images resulted from the work of this thesis, but also new ideas have surfaced. The final chapter of this thesis therefore contains a short discussion on some directions for future research. These include possible adaptations of the proposed techniques as well as 19 other directions of investigation that now seem worthwhile pursuing. 20 Chapter 2 A Study of the Phase Wrap Artifact in MRI: A New Method of Phase Unwrapping 2.1 Understanding Phase Maps in MRI As described in the previous chapter, the MRI signal is composed of real and imaginary components or equally, it can be described in terms of a magnitude and phase. The phase of the MRI signal can be displayed as what will be referred to as a phase map. Important information, such as field inhomogeneity or the velocity of blood flow, can be encoded in the phase of the MRI signal [18]. Chapter 4 will have more details on velocity phase encoded MRI. Here, the focus is not on how such an image is obtained, but rather how to extract the useful information given such an image. Although phase maps can contain very important information in MRI and many other fields of study, they are not well exploited. This is due to the tricky nature of the phase. Any cyclic waveform can be described at any instant, in terms of its magnitude and phase. However, the phase information only has relative significance since it indicates the stage the waveform is at relative to some reference. In MRI, this ambiguity has made it so that the magnitude of the complex MRI signal is the information that is primarily used for clinical diagnostics. This will inevitably result in the loss of important information which is and can only be encoded in the phase of the signal. In this chapter, MRI phase maps are studied. In particular, the most common artifact associated with these maps and greatly limiting their use, called phase wrapping, is explored. The characterisation of such phase maps as well as a better appreciation of the difficulty presented by the phase wrap artifact in MRI, led to a new method of phase unwrapping. This chapter begins by describing the characteristics and properties of wrapped phase 21 maps in general and then, MRI phase maps are discussed. Such maps must be well under-stood in order to appreciate the problem of phase unwrapping. In section 2.2, the problem of phase unwrapping is presented along with a review of some phase unwrapping methods that are presently available. In section 2.3, the motivation that led to the new unwrapping method is given along with a description of the steps involved and some results. Further-more, the method is tested for its robustness to noise and an extension of the work is proposed in section 2.4. 2.1.1 The Nature of Phase Information One way to view the phase of any complex signal, is that it represents a rotation, with direction and amplitude. However, given any complex data, the phase can only be derived as modulo 2TT which gives a scalar value (known as the principal value) contained in a given range, usually [—7r,7r) and the true phase rotation cannot be known unambiguously from a single scalar value. For example, rotations of 7r/2 and —ZTX/2 are indistinguishable since —37r/2 falls outside the unambiguous range [—7r,7r) and is hence wrapped around to give a principal value of T T / 2 . This wrapping procedure consists of adding or subtracting a multiple of 2TV such that the result is contained within [—7r,7r). If phase maps are to convey any useful information, the exact direction and amplitude of rotation of the signal in one pixel with respect to that in another must be known. Any error made in assigning a phase value to a pixel will result in an erroneous quantification of the underlying physics. For example, flow may be misinterpreted as going in the wrong direction and with the wrong magnitude (off by a factor of three) if a rotation of — 37r/2 is not distinguished from a rotation of -\-TT/2. When phase maps contain no wrappings, each principal value indicates the exact am-plitude and direction of rotation. These types of maps will not be further discussed since they already contain readily accessible information, just like magnitude maps. 2.1.2 Phase Loops : Poles Wrapped phase maps may contain two types of borderlines called fringelines and outlines [19]. It is important to distinguish between these structures when trying to understand the information contained in a phase map. Fringelines are borderlines between two adjacent pixels where phase wrapping appears to occur since it denotes a +7r/ — TI boundary. However, the principal phase values on either side of a fringeline may in fact be the true phase values and hence represent a physical reality. In such a case, this fringeline will not have resulted from a wrapping procedure and will also be a cutline. Outlines are borderlines between two adjacent pixels where the signal has undergone a relative rotation of more than TX (in either 22 direction). Outlines are often more difficult to discern than fringelines since they do not necessarily follow an obviously contrasting line. In most cases, fringelines and cutlines will overlap in some areas but not in others. When no cutline is present, it is straightforward to determine the exact value of the phase at any point by adding or subtracting a multiple of 2ir each time a fringeline is crossed. This technique is known as fringe counting or fringe scanning [20]. This is due to the fact that a fringeline in such a case has been caused by the wrapping procedure so the true phase can be determined from the principal phase value by an analogous unwrapping procedure. This is only possible when no cutlines are present so the true phase varies by less than TX between adjacent pixels. Such cases occur when a fringeline is closed with no open ends as in Fig.2.1. However, the difficulty in phase map interpretation arises when a fringeline has Figure 2.1: (a) Phase map simulation with closed fringeline and no poles, (b) The profile of (a) along the dotted line. The wrapped signal can be unwrapped in this case by adding 2TX when crossing the fringeline from bright to dark. The true (unwrapped) phase difference between points A and B is less than ir. two open ends due to the presence of a cutline as in Fig.2.2. In these cases, the four pixels forming each of the ends result in phase loops known as inconsistencies or residues [11]. The term pole will be used here to emphasize that the values +2TX and — 2TX associated with such four-pixel structures can be of either sign. In fact, one end of an open fringeline is a positive pole and the other is a negative pole. Such a pair will be referred to as a dipole structure (see Fig.2.2). If the boundary conditions are such that the phase map is connected at both ends, the image has the topology of a torus as shown in Fig.2.3. This is the realistic case for 2DFT MRI, described in the previous chapter, since the complex image is the result of a 2D Fourier transform of the 2D complex NMR signal. In such a case, there are always an equal number / c l o s e d fringelin (a) (b) 23 (a) (b) Figure 2.2: (a) Phase map simulation with open fringeline and a dipole structure, (b) The profile of (a) along the dotted line. The wrapped signal cannot be easily unwrapped because errors will be made at the poles. The cutline causing the poles is the borderline between grey and dark region, assumed to terminate at the poles on either end. This cutline will cause the unwrapping to be path dependent. The true phase difference between points A and B is greater than n. of positive and negative poles since all open fringelines have both ends defined within the phase map. For all further discussions, all poles considered belong to a dipole (no monopoles [21] are present). An interesting feature of the dipole structures is that the pairing up of positive and negative poles can be affected by a global phase offset. Although a global phase offset will not affect the pole positions, it will affect their connections as the fringelines are shifted. In other words, a positive pole may no longer connect to the same negative pole for a given phase offset. This is demonstrated in Fig.2.4. Given a two dimensional phase map (f)(x,y), a pole, p , can be defined mathematically as a 2 x 2 pixel loop for which the curl of the two dimensional wrapped phase gradient G^ is non-zero: P = V x G ^ 0 (2.1) where : G 0 = W[&4>x{x,y)] i + W[A<f>y{x,y)] j (2.2) with A(f)x(x,y) and A(f)y(x,y) defined by: A(f)x(x,y) = <j)(x + l,y) - <j>(x,y) A(j)y(x:y) = cj)(x,y + 1) - cj)(x,y) (2.3) and where W[-] is the wrapping operator which gives a result contained within the range [—7T, IT) by adding or subtracting a multiple of 2ir to whatever is in the square brackets. 24 no boundary conditions rectangle one boundary condition A j k C cylinder two boundary conditions torus Figure 2.3: The boundary conditions in MRI are such that the phase map is connected at both ends having the topology of a torus. More specifically, p is given by: p = {-W[A(j>v(x,y)]-W[A(l>x(x,y + l)] +W[A<j>y(x + 1,y)] + W[A(j)x(x, y)}} k. (2.4) To evaluate this, the following relations of the wrapping operator can be used: -W[Ac/>ab} = W[-Act)ab\ (2.5) W{A<j>ab] + W[A(fibc} + W[A(}>cd] = W[A4>ad] + 27rm (2.6) where : f -1 for W[Acf>ab} + W[A<f>bc] + W[A<f>cd\ < - T T ?n=\ +1 for W[A4>ab\ + W[Afoc] + W[A4>cd] > TT (2.7) [ 0 for \W[A(f)ab} + W[A(t>bc) + W[A(t>cd]\ < vr and where the notation: Acj)ab indicates </>a — cj)b. The expression for p thus yields: p = {W[-A<f>y(x,y)] + W[-A<j>x{x,y+1)} +W[A<f>y{x + 1,7/)] + W[A<f>x(x, V)}} k = {W[-A<f>x(x, y)] + 27rm + ^ [ A ^ ^ , y)]} k = {-WiA^ix, y)] + 2nm + W[A<j>x{x, y)}} k = 2irmk (2.8) where again m = ±1,0. Eq.2.8 and Eq.2.1 give the result : p = ±27r k , corresponding to positive and negative poles. Hence, any 2 x 2 pixel loop will be a pole if m = ±1 and it will 25 ( a ) (b) ( c ) Figure 2.4: Phase maps derived from a single random complex field with various phase offsets: (a) 0°, (b) 45° and (c) 90°. The opposing poles of a given dipole, o, change connection from (a) to (b) whereas those of another dipole, • , change connection from (b) to (c). not be a pole only for the case m = 0 which, according to Eq.2.7, results only if: This equation can therefore be referred to as the no pole condition. This formulation is analogous to that of conservative and non-conservative fields in mechanics. Mathematically, a conservative field F is one which can be expressed as the standard gradient of a scalar function V: F = \jV. This results in a path independent integration : / F • dr, where dr is a differential spatial vector along the integration path in three dimensions. Equally, this results in a vanishing curl : V x F = 0 since the curl of a gradient field is always zero. A standard gradient field derived from an MRI phase map should therefore be a conser-vative field where the properties of path independent integration and vanishing curl apply. However, the gradient field is defined as the wrapped gradient field of a scalar phase map cf)(x,y). It is clear from Eq.2.1, that the wrapping operation may produce regions in G ^ where poles are present (a non-vanishing curl) and henceforth the path independent integration property is no longer valid in these regions. The gradient field G ^ can thus be considered a non-conservative field if poles are present. Furthermore, for any vector field F , the following relationship holds: where again dr is the differential spatial vector in three dimensions, and dS is a differential surface vector perpendicular to S, the surface enclosed by the line integral, and indicates \W[-A</>y(x, y)] + W[-A<f>x(x, y + 1)] +W[A(/>y(x + l,y)\ < TT (2.9) (2.10) 26 a surface integration. For a two-dimensional field such as G^(x,y), the curl is in the z direction as well as dS since S is in the (x, y) plane. Using the definition of the pole given in Eq.2.1, Eq.2.10 gives: where pz is the z component (only component) of the pole and dAz indicates integration over the area of the (x,y) phase map enclosed by the line integration. Eq.2.11 implies that the integration of G^,(x,y) along a closed path will be equal to the sum of the pole values, pz, contained within that region. If that integration is zero, either there are no poles within the region or there are an equal number of positive and negative poles. The converse is true as well. Therefore, closed-path line integration of the phase gradient map is path independent as long as no poles or an equal number of positive and negative poles are enclosed by the path. In such a case, the integration is considered consistent since any path will yield the same result. The cutline phase unwrapping methods ensure consistency in the integration of phase maps by placing cuts between the dipoles, guiding the integration, as will be.discussed in section 2.2. As defined in the previous section, poles are present at the ends of the open-ended fringelines. Since poles result from a true phase difference of more than 7r between the signal in two adjacent pixels, the borderline between these two pixels is a cutline. In fact, the shortest possible cutline is a single pixel edge in length with positive and negative poles at either end sharing two common pixels, one on either side of the edge. Therefore, open-ended fringelines indicate the presence of a cutline which can be assumed to begin and end at the poles but the rest of it is not easily detectable [19]. Open-ended fringelines or cutlines may be caused by at least three different factors in MRI and hence there are at least three different sources of poles. One such source is noise. Because the noise in MRI is a small random vector added to the complex signal, the noise will determine the phase in areas of low relevant signal. Since this type of noise is random, the phase in these regions will randomly fluctuate and an open-ended + 7 r / — IT boundary may result. Fig.2.4 shows a phase map corresponding to a smoothed noise field. Several poles are indicated by the circles and squares. In addition, Fig.2.5 shows a zoomed section of a purely random phase map where several poles are circled and their corresponding fringelines are indicated by the grey dashes. In fact, it can be shown statistically that, for a purely random phase map 4>random{x, y), with torus topology and N x N pixels, there will be approximately (N x N)/3 poles in all. (2.11) 2.1.3 Sources of Poles in MRI 27 Figure 2.5: Zoom of a random phase map showing several poles (circled) and their fringeline connections (dashed borderlines). Furthermore, since such boundary conditions will lead to an equal number of positive and negative poles, there will be approximately (JVx/V)/6 positive poles and (NxN)/6 negative poles. To prove these statements, let Pw[A<j>] be the probability distribution for the wrapped phase difference of any two four-connected pixels (connected to the neighbours that share a common edge). Since the difference is wrapped, any value for — TT < W[A<f>] < 7r will be as likely to occur so Pw[A<f>] will be a step function which, when normalized, can be written as: Using the no pole condition expressed in Eq.2.9, let W[A^i], VP]A02] and VF[A03] be the three wrapped differences involved in Eq.2.9. These differences can be any three of the possi-ble four associated with a 2 x 2 pixel loop. The probability of being a pole, either positive or negative, P±, will thus be equal to the probability of the sum of these three phase differences being equal to a value outside the range (—7r, 7r). Here, the unambiguous range is considered to exclude the —7r end-point for simplicity. In other words, given the normalized probability distribution for the sum of three wrapped phase differences, P{w[A4>l\+w[A4>2}+w[A(t>3}) written as P(Wi+w2+w3) f ° r simplicity, P± can be calculated as follows: where d{W\ + W2 + W3) indicates integration over the variable W\ + W2 + W3 which is simply the sum: W[A0i] + W[A0 2] + W[A0 3]. Since the wrapping operator makes it such that W^A^j], iy[A</>2] and W[A0 3] are all independent of each other, P(Wl+w2+w3) can be easily calculated given Pw[A<j>{\ = Pw[A<t>2] = ± for - T T < W[A<j>] < 7T f/ for W[A(j)} < - T T k W[A</>] > TT (2.12) P± = (2.13) 28 Pw[A<p3] = Pw[/\<t>] by simple convolution (*): P(w1+w2+w3) = Pw[A<j>] * P\v[A<t>] * Pw[A<t>] (2-14) where PW[A<J>] is given in Eq.2.12r The result of the convolution, after normalization, gives: P(E^ = \ ^ ^ f o r 7 r < | | E ^ | | < 3 7 r (2-15) [ 0 for \\T,W\\ > 3TT where £ W = W1 + W2 + W3 = W [ A & ] + T^[A^2] + W[A(/)3}. A graph of this function is shown in Fig.2.6. Figure 2.6: Probability plot for the value of the sum of three wrapped phase differences, P^w)- ^ n e function is given by Eq.2.15 and can be used with Eq.2.13 to obtain the result: P± = 1/3 as indicated by the shaded area. This result indicates that given a random phase map made up of N x N pixels, approximately (N x A r)/3 poles will be present if N is large. Evaluating Eq.2.13 using Eq.2.15 gives the result P± = 1/3 represented by the shaded area of the curve in Fig.2.6. It is clear from the symmetry of the graph shown in Fig.2.6 that if P+ and P_ are the probabilities of being a positive pole and a negative pole respectively, P+ = P- = P±/2 — 1/6. If a random phase map, with torus topology, is made up of N x N pixels, there will be a total of iV x N potential poles hence P± predicts that approximately (JV x jV)/3 poles will be present if N is large. This result has been verified for random phase maps. 29 In a purely random noise phase field, poles randomly appear and it does not make sense to assume that the poles are due to wrapping, so the principal phase value can be considered the true phase value. Applying this rationale to MR images, in noisy regions, all fringelines can be considered cutlines and no unwrapping needs to be done. Another source of poles is spatial undersampling of the phase. This will result when the phase sensitivity and spatial resolution of the imaging sequence are such that two adjacent pixels have a true signal rotation of more than TX. In theory, most physical fields mapped by the phase of the MRI signal are continuous. This means that, in theory, a fine enough spatial sampling will always eliminate such large phase differences between adjacent pixels. However, in practice this is not the case since many other factors will limit the spatial resolution of the MR image as mentioned in Chapter 1. Fig.2.7 demonstrates how poles appear when the sampling resolution of a continuous function is made very coarse. ( a ) <fc>) ( c ) ( d ) Figure 2.7: Four phase maps are shown with corresponding central horizontal profiles under them. A continuous function in (a) is sampled with different resolutions: (b)16 pixels, (c) 32 pixels, (d) 42 pixels. Poles only appear in (d) as indicated by the circles and the larger than TX phase jumps (|A0| > TX) between adjacent pixels. There are at least two different cases in MRI, when these types of undersampling poles can be seen. The first case is that of velocity encoded imaging, where the flow is encoded in the phase of the signal. This is the central theme of Chapter 4 where more details are given as to how such images are produced. An example of such a case is shown in Fig.2.8 and others can be seen in Chapter 4. Such poles are usually avoided by reducing the flow sensitivity of the imaging sequence. However, the reduced flow sensitivity can cause the loss of other important information 30 Figure 2.8: Velocity encoded phase map of blood circulating in the heart. The fringeline is open-ended. One pole of either sign is present at the fringeline ends. by reducing the range of blood flow velocities that can be distinguished. Also, it is not always possible to predict a priori the expected maximum blood flow, specially for diseased patients, hence undersampling poles could still result. If a phase unwrapping method could be used to overcome this limitation on the flow sensitivity, velocity encoded MRI, and, as a consequence, the dynamic flow imaging of Chapter 4, would become more practical and useful. The other case of undersampling poles is due to main field inhomogeneity. This inho-mogeneity may result from magnetic susceptibility differences such as those that occur at tissue-air interfaces. An imaging sequence that uses a gradient echo is sensitive to this inho-mogeneity and thus may result in a discontinuous phase map, where the phase of the signal in adjacent pixels varies by more than ix. Very nice gradient echo images demonstrating such open-ended fringelines with undersampling poles were obtained by Reichenbach et al. [22] (see figure 5 of their paper). The main distinction between poles arising from undersampling and those arising from noise is in the length of the fringelines connecting the dipoles. In regions of low signal-to-noise ratio (SNR), dominated by noise, the phase coherence is also low since the values of the pixels are random. This means that a -T-IX/ — ix boundary appearing somewhere does not affect the probability that adjacent pixels will also have a +ix/ — ix boundary. The result is that in such regions, a fringeline will not tend to be very long. In contrast to this, in regions of high SNR, the phase coherence is high since the phase is dominated by the relevant signal which is in general, a slowly varying continuous field. In such a case, if there is a +1x1 — ix boundary at one location, there is a higher chance that the regions around it will also contain a +ixj — ix boundary. The result is that in regions of high SNR, the fringelines tend to be longer than those in the noisier regions. 31 2.2 The Phase Unwrapping Problem 2.2.1 Review of the Phase Unwrapping Problem In order to access the information contained in a wrapped, principal value phase map, it is often necessary to first unwrap it successfully. The problem of phase unwrapping can be stated mathematically as finding fame given ^wrapped a n d the fact that 4 > w r a p p e d — <f>true + 27rn = W[0true] where W is the wrapping operator, n is an integer and —TT < 4>wrapped < T T - Therefore, the problem is basically one of finding the correct value for n at every pixel. If the true phase difference between adjacent pixels is contained within the unambiguous range [—TT,TT), this problem becomes quite trivial : a line integration of the gradient (phase difference between adjacent pixels) along any path, will result in an unwrapped phase map. Furthermore, fame at every pixel can be determined from the value of the resulting unwrapped phase map if n is known for the starting point of the path of integration. For example, if the starting point for the gradient integration is taken to be a pixel for which n = 0 and (j>wrapped = <f>true, no correction will be necessary and the resulting unwrapped phase map from the integration represents the true phase map. If n ^  0, all the pixels in the resulting unwrapped map will be shifted by 2-nn where n is the number of cycles by which the initial pixel value is wrapped. This is the most basic and straightforward method to phase unwrap and is what is commonly used to solve this problem. However, if the phase of adjacent pixels does not obey the — TT < Acf)true < TT condition, such a straight forward solution will fail. This is when the problem of phase unwrapping becomes challenging and the focus of many publications. The reason the —7r < A(j)true < 7r condition must be met is due to the fact that this is implicitly assumed when integrating the gradient since the gradient can only be known modulo 2TT. When this condition is not met, the result of the integration becomes path dependent. In such cases, integrating across neighbouring pixels for which A(j)true is outside the range [—IT, TT) will cause errors to propagate across the rest of the unwrapped image which will then be off by some unknown multiple of 2TT. This is what sometimes appears as "streaks" if the integration is performed along rows or columns. Fig.2.9 shows a typical phase map image that results when the straightforward unwrapping approach fails. Most of the work in phase unwrapping has been done in other fields such as optical engineering. In particular, most new phase unwrapping developments are for synthetic aperture radar (SAR) interferometry. However, the applications can easily be extended to MRI. Good reviews of phase unwrapping methods are given by Judge and Bryanston-Cross [23] and Ghiglia and Pritt [11]. 32 Figure 2.9: (a) Phase map representing flow in a longitudinal direction (b) Phase map with streaks resulting from conventional phase unwrapping algorithm 2.2.2 Review of Phase Unwrapping Methods The present methods of phase unwrapping can be grouped into two major categories: (i) minimum normalisation methods and (ii) path-following methods [11]. In this thesis, the new method of phase unwrapping proposed belongs to the second category of phase un-wrapping methods, the path-following methods. The goal of this work is not to compare its performance to that of other methods although the motivation for the proposed method is based on the limitations of other methods. Comparisons of the performance of different methods can be found elsewhere [11, 24]. The following is a general review of phase un-wrapping methods. Only the most important aspects of relevant phase unwrapping methods is given in order to motivate the development of this new method while setting its place amongst the long history of phase unwrapping work. The recent textbook written by Ghiglia and Pritt [11] provides a more thorough review of the present status of phase unwrapping with the details, algorithms and even software of the many different phase unwrapping techniques available. The first category of phase unwrapping methods, the minimum normalisation methods, impose constraints on the solution. They do this in a mathematically formal way where, in general, some sort of fitting procedure is performed to best match the local derivatives of the unwrapped phase map solution, to the measured derivatives in the wrapped phase map. The true phase map is therefore derived by solving a least-squares problem. In this category are certain surface-fitting methods and several Poisson equation based methods [11]. The problem with these methods is that they do not guarantee a streak-free phase map since they do not explicitly attempt to block the propagation of phase errors. Furthermore, an image such as that shown in Fig.2.8 with a large discontinuity could not be well modeled by such methods. Iterative variations of these methods have produced better results but are limited by their inefficiencies. Also, in many cases, these methods are too computational to 33 be practical. Since the new method for phase unwrapping, proposed in this thesis, follows such a different approach from these minimium normalisation methods, they will not be further discussed here. For more details on these methods of phase unwrapping, see [11]. The second category of phase unwrapping methods, the path-following methods, deal with the idea of guiding the integration of the gradient map so as to avoid any error prop-agation. More specifically, these methods rely on integrating the gradient map along a particular path which avoids using any true phase differences which are outside the unam-biguous range [—TT, TT) to unwrap. In this category are the disk growing method [25], cut-line methods [26, 27, 28] and region grow methods such as minimum and maximum spanning tree methods [29, 30, 31]. The temporal phase unwrapping proposed by Xiang [32] can also be included in this category since it proposes to integrate the phase differences along another dimension in order to avoid integrating across a true spatial phase jump outside the range [—TT, TT). This method will only work if the familiar condition : — TT < Acfi < TT is obeyed for pixels in temporally adjacent frames. This condition along with motion between temporal frames will further restrict and complicate the application of this procedure [32]. Path-following phase unwrapping methods often result in streak-free phase maps. This is because their main goal is to block error propagation although the result may not always be correct. As long as the poles are connected or blocked somehow from entering the unwrapping process, any integration path, starting at the same initial pixel, will give the same result or be considered consistent. This result will however be dependent on the exact placement of cuts or criterion used to grow. As appreciated by Ghiglia and Pritt [11], the choice of cuts is not obvious, some criterion must be used. Given the description of fringelines and cutlines in section 2.1, it is clear that in order to obtain the true phase map, unwrapping should be blocked along the cutlines. This is because it is along these borderlines that true phase differences are outside the unambiguous range. This reveals the choice of terminology since these cuts should be used to block the unwrapping along the cutlines. However, most cut-line phase unwrapping methods connect the poles and block the unwrapping but they do not explicitly look for cutlines. Instead, they use arbitrary criteria like minimizing the cut lengths [11, 23] to determine the placement of the cuts. The reason so many of these methods seem successful is only because they deal with noise induced poles. As discussed previously, noise poles will have short fringelines/cutlines associated with them so the minimization of length criterion is justified. As long as all the poles in a phase map are the result of noise, cuts connecting these along shortest paths result in good representations of the true phase map. This is due to the fact that even if the cuts do not follow the cutlines closely, only small areas of the phase map are affected. However, most of these methods will fail to result in true phase maps if there are long range dipoles 34 due to undersampling by either phase encoded velocity or magnetic field inhomogeneity. It is in these cases, with longer cutlines, that any error in the placement of the cuts will affect a larger region of the map. No phase unwrapping method which explicitly places cuts has been able to deal with such dipoles very successfully. Region grow methods, which guide the unwrapping without explicitly determining the placement of the cuts, have been shown to successfully unwrap phase maps with large order open-ended fringelines [28]. This is to be expected if the criterion used to guide the growing, such as minimum phase gradient, results in circumvention of the cutlines. In these cases, there will be effective cuts implicitly placed along the cutlines since there is a direct correspondence between the region grow path and the cuts [30]. The limitation of these methods is the speed [23]. The correspondence by An et al. [33] introduces a new implementation of a minimum spanning tree phase unwrapping method where the speed has been greatly improved, reducing it from order 0(n 2) to 0(nlog2n). However, these methods are also difficult to track since the effective cuts are not as easy to see. This may make the result more difficult to understand. Ghiglia and Pritt [11] present a phase unwrapping algorithm called : mask cut algorithm. This algorithm was developed as a hybrid method combining the ideas of cut-line methods with region grow methods as was done by Flynn [34]. Region grow methods use some other information (quality map) to help guide the growing. Cut-line methods place explicit cuts between the dipoles to avoid path-dependent results. The idea behind the mask cut algorithm is to choose the cut placements according to a quality map so that not only the poles are blocked from the unwrapping path, but additional information is also incorporated. A similar motivation led to the cutline method presented in this thesis. 2.3 The New Cutline Phase Unwrapping Method 2.3.1 Method The understanding of the source of the poles led to a new idea for finding the cutlines ex-plicitly. The idea relies on a score map generated by repetitive tracking of the fringelines, connecting the dipoles, for several phase offsets. This score map is then eroded (thinned) so as to obtain the cuts. Kramer and Loffeld [19] suggested a similar (perhaps more compli-cated) fringline tracking procedure for cutline detection. However, in contrast to the method presented here where the fringeline tracking is used to produce a score map, Kramer and Loffeld used fringeline tracking to directly extract portions of cutlines. As pointed out in their discussion, this led to several limiting factors for their method. Also, their paper con-centrates on cutline detection without attempting to phase unwrap which can be a problem 35 if the entire outline is not found as was the case for their results. A new method of finding the location for the cuts can be derived from the understanding of the phase map as described in previous sections. This new phase unwrapping method will be referred to as a pole-guided-cutline phase unwrapping method. As described above, poles are point-like (2 x 2) inconsistencies present in phase maps. Nonetheless, it is well appreciated that removal of the poles (i.e. replacement of the pole pixels by a fixed value) alone is not sufficient to guarantee a streak-free phase unwrapping result [35]. In fact, removal consecutively creates adjacent poles along a path connecting positive and negative poles until they superimpose and cancel in a process called annihilation. This phenomenon is the basis for the cutline detection that is proposed here. The first step is to distinguish between noise dipoles and those which are longer-range and may be due to undersampling. However, this is not a necessary step since the noise poles could be treated like the longer-range poles without affecting the result. The distinction between these two types of poles is simply used to make the method more efficient since it is not necessary to spend as much time trying to determine the cutlines associated with noise dipoles. Instead, the fringelines can simply be assumed to be cutlines. The distinction between the types of poles is made by tracking the fringelines and then classifying the poles according to their relative fringeline lengths. Fringeline tracking is done by consecutively chasing the poles. At each step a pole is removed (all four pixels are replaced by zero phase value) thus creating a new one in one of the eight-possible adjacent 2 x 2 loops containing at least one of the original pole pixels. This process continues until no new poles are found at which point dipole annihilation has occurred. The result of the tracking will be a two-pixel-wide line of phase equal to zero which delineates the fringeline. A zoomed version of the phase map can be used for these steps since otherwise, removal of a pole displaces all surrounding poles. If all fringelines in the image are tracked simultaneously, the result will be an initial reduction in the total number of poles as the noise dipoles are annihilating. When the total number of poles no longer decreases with pole removal, the leftover poles can be classified as longer-range and their corresponding fringelines treated accordingly. This will allow for simple noise extraction without any thresholding. Furthermore, if this method results in any noise poles being mistaken for long-range poles and treated accordingly this will not affect the result. It will only increase the amount of unnecessary calculations. The cuts connecting noise poles can be determined from a single fringeline tracking procedure since the fringelines are cutlines in the noise. However, determination of the cuts connecting the other dipoles will require a lot more work. The longer cutlines for the other long-range dipoles are determined by repeating the fringeline tracking for several global 36 phase offsets and superimposing the resulting lines. The idea, as also stated by Kramer and Loffeld [19], is that those parts of the fringeline connecting a given dipole that do not vary much with phase offset are part of the cutline. This is due to the fact that the smaller the phase gradient in a region, the more the phase offset will cause the fringeline to vary in that region. The converse is true as well. The sought cutline will then be visible in the superimposed map, which can now serve as a score map, as the most highlighted connection of the dipoles. This line can be extracted by first region growing from one pole along the regions of higher signal (more overlap) until a pole of opposite polarity is reached and an island is formed. This island may have several "holes" which are in turn filled up, row by row. Then, an erosion procedure described by Xiang et al. [36]. is used to convert the island into a thin cut. This procedure is called the score-guided erosion since it uses the information from the score map in order to determine in which order the pixels will be eroded (removed). It begins by looking at all the pixels that make up the edges of the island. The pixels with lowest corresponding score values (according to the score map) are sequentially eroded as long as they do not affect the connectivity of the remaining island. At each step, new edge pixels are included as possible candidates for erosion. The connectivity check insures that the final result will be a line connecting the two poles of opposite polarity. Alternatively, the mask-thinning algorithm described by Ghiglia and Pritt [11] or Flynn [34] for their mask cut phase unwrapping algorithms could be used. A more efficient implementation of the described fringeline tracking and phase offsetting is to change the phase of the value used for pole removal with no need to do a global phase offset. In other words, instead of setting the pixels to zero phase, they are consecutively set to another value n • 9STEP where n = 0,1.. .N is the number of overlapping fringelines and N = 27r/estep. Once the cuts are all determined, unwrapping can proceed in any fashion as long as it is blocked at the cuts. The implementation used in this paper is to unwrap into but not past the cuts. If some pixels are surrounded by cuts or if the cuts are thick (greater than three pixels wide) in some areas, there will be some left-over, unwrapped pixels at the end of all the unwrapping. This will occur mostly in noisy regions. These pixels can be unwrapped at the very end, layer by layer. The result will be some unwrapping in the noise but it will not corrupt the result. Flynn's phase unwrapping guided by a quality map [34] is very similar to the phase unwrapping presented here. They differ in that the quality map used by Flynn is magnitude-based whereas the score map proposed here is only phase dependent. The score map is the core contribution of this work. It emphasizes high phase gradients only in regions which are relevant for dipole connections. The pole-guided-cutline method can be seen as 37 a generalization of Flynn's method to be used whenever a magnitude-based quality map is not readily available. 2.3.2 Results In this section, some results from the previously described phase unwrapping method are presented. Fig.2.10 shows four original phase maps with the successful results of the pole-guided-cutline phase unwrapping method. It is interesting to note that the cutline method (a) (b) Figure 2.10: Four velocity encoded images I, II, III and IV. The original phase maps are shown in (a). Unwrapped results from the pole-guided-cutline method are shown in (b). works even when there are no long-range dipoles (as in Fig.2.10 II , where all the high signal fringelines are closed). Such a case is very simple to unwrap since no score map and long-range cutline determination is necessary. The fringelines in the noise can simply be used as cuts for the unwrapping. Fig.2.11 shows the different stages of the cutline method for one particular example (Fig.2.10 I). The original velocity encoded phase map and the resulting unwrapped map are shown in (a) and (b) of Fig.2.10 I and Fig.2.11, respectively. Once the classification of poles has been done, the location of long-range dipoles can be used to zoom in on the relevant part of the original phase image. This relevant region is shown in all images of Fig.2.11. In (a), arrows point to the long-range dipoles with white indicating positive and black negative. Many of these dipoles are clearly undersampling dipoles since they occur at the ends of 38 Figure 2.11: (a) Velocity encoded phase map with long-range positive (white) and negative (black) dipole locations indicated by arrows, (b) Unwrapped result using the pole-guided-cutline method on (a), (c) Noise dipole cuts, (d) The score map obtained by overlaying long-range dipole fringelines for various phase offsets, (e) Cuts resulting from the score-guided-erosion (SGE) performed on (d) given the dipole locations indicated in (a). high signal (flow) fringelines. However, a fringeline in the noisy region (right) was seen to have a similar length and hence the associated dipoles are considered in the long-range cutline treatment. As previously mentioned, this will not affect the result but only extend the processing time. Noise dipole fringelines (cuts) resulting from fringeline tracking are shown in (c). The interesting score map that is obtained by overlaying the dipole fringelines corresponding to various phase offsets is shown in (d). This map is used to score the pixels when performing the region grow and score-guided-erosion (SGE) described by Xiang et al. [36] since the brightest pixels should be used to find the cuts. Image (e) gives the resulting cuts which can be seen to follow the brighter parts of (d) while connecting the dipoles pointed out in (a). The cuts shown in (c) and (e) are used to obtain the unwrapped result in (b). 2.3.3 Study of Noise Effects This pole-guided-cutline phase unwrapping method relies on the score map to obtain the cuts but even if there is a lot of noise corrupting the image, the score map produced should 39 still be useful since all the dipoles are connected at every phase offset step. The noise is not expected to disrupt the final connectivity between dipoles as was the case for [19]. This is also in contrast to the mask cut method [11, 34] which relies on a good quality map and fails when the quality map is too corrupted. To test the robustness of the pole-guided-cutline phase unwrapping method presented in this thesis, the effects of noise were studied. A simulated phase map was produced with varying levels of Gaussian noise added to the corresponding real and imaginary maps. The signal-to-noise ratio (SNR) is calculated, in a flat region of the magnitude image, as the average signal divided by the standard deviation. The unwrapped results are shown in Figs.2.12 and 2.13. SNR = 17.5 SNR = 6.0 SNR = 3.3 SNR = 2.9 Figure 2.12: Simulations of a wrapped phase map with no cutlines is shown in row (a). From left to right, the signal-to-noise ratio (SNR) deteriorates as indicated in the top row. Row (b) shows the result from applying the pole-guided-cutline phase unwrapping method to the corresponding image in (a). For the first simulation shown in Fig.2.12, there is no cutline so all poles are due to noise. The pole-guided-cutline method result was not affected by the noise. In Fig.2.13, there is a long-range cutline present along the bottom of the circle and again, the pole-guided-cutline method result was not affected by the noise. In conclusion, the detection of the long-range cutline connecting undersampling poles is immune to the effects of noise, making the pole-guided-cutline method very robust. 2.4 Three-Dimensional Phase Unwrapping 2.4.1 Limitations of Two-Dimensional Phase Unwrapping An important limitation of any two dimensional phase unwrapping method is the inability to unwrap a fully wrapped region, with no detectable fringelines. If the three dimensional 40 S N R = 1 7 . 5 S N R = 6 . 0 S N R = 3 . 3 S N R = 2 . 7 Figure 2.13: Simulations of a wrapped phase map with a cutline along the bottom is shown in row (a). From left to right, the signal-to-noise ratio (SNR) deteriorates as indicated in the top row. Row (b) shows the result from applying the pole-guided-cutline phase unwrapping method to the corresponding image in (a)-information is available, either for cine velocity encoded images [37] with (x, y, t) infor-mation or for field inhomogeneity mapping with (x, y, z) dimensions, unwrapping in three dimensions is expected to improve the results. One example of such a fully wrapped region in the x — y plane is that of a fully wrapped vessel in cine velocity encoded MRI (see Chapter 4) as seen in the lower part of Fig.2.10 III. For such a case, no poles will be detected in the x — y plane so there will be no fringeline tracking in the signal. This is the same as what happens in the closed fringeline case (such as in Fig.2.1) which can be successfully unwrapped by any straightforward method if the noise is treated somehow. However, for the fully wrapped vessel there is no good (unwrapped) signal to provide the information needed to successfully unwrap. If the fully wrapped vessel is treated in three dimensions, the problem can be solved. The unwrapped signal from temporal frames adjacent to a fully wrapped vessel frame can provide the necessary information for temporal phase unwrapping [32]. This direction for the unwrapping will be favoured when avoiding large phase jumps in any of the three dimensions. This is because although no poles are found in the x — y plane, a ring of poles and hence large phase jumps in the (x, t) and (y, t) dimensions can be detected all along the vessel edge as shown in Fig.2.14. This ring of poles will lie parallel to the x — y plane, at the border between ti and t2 as well as between t2 and £3. These poles will block any spatial unwrapping from the noise region into the vessel while allowing for temporal unwrapping within the vessel. The pole-guided-cutline unwrapping method avoids areas of large phase jumps when 41 V p r o f i l e s H2 p o l e s p o l e s in ( x l V ) \niyJX 7T -TC B TC • r e 71 -TC (a) , A R < Tt B wy-w|<* "7 (b) Figure 2.14: (a) Sketch representing the phase map of a wrapped vessel at three different moments in time: t\, ti and t$ with three dimensional poles shown by the dashed rings. The vessel is fully wrapped at t-2- (b) The profiles shown in the first column represent the phase (p as a function of y. The second column shows there will be no poles detected in any of the three (x,y) phase maps. The third column shows that there will be poles present in the {y,t) map at the edge of the vessel. These poles will arise due to the large phase jumps between the different time steps for the same point in space (as shown for point A). If all possible poles are searched for in three dimensions, the result will be a ring of poles lying along the edge of the vessel, parallel to the {x,y) plane, at the border between ti and t-2 as well as between £2 and £3 as shown by the dashed lines in (a). performing the phase unwrapping integration. Extended to three dimensions, it should favour the time dimension as a path to enter the fully wrapped vessel since it is spatially surrounded by large phase jumps. The formulation of section II can easily be extended to three dimensions since the curl in Eq.2.1 can be applied to a three dimensional field G^,(x,y,t). The result will be a pole vector with discretized components (either 0, +27T or —2TT) in each of the three directions. Three dimensional poles can hence be searched for in a three dimensional phase map where fringelines are now fringesMr/aces and cutlines are cutsurfaces. In theory, the same steps as those used for the two dimensional cutline phase unwrapping method can be used to find the surfaces which should block the unwrapping. The third dimension also has the potential of providing additional information about the dipoles. In two dimensions it was shown that fringelines will change dipole connections 42 depending on the phase offset (see Fig.2.4). For the pole-guided-cutline method, fringelines in the noise are treated as cutlines hence the noise dipoles are connected according to the fringeline connections with no phase offset. For the longer range, high signal dipoles, correct connections between the dipoles and placement of the cuts relies on the presence of a larger phase gradient that is outlined by the score map. In both cases, assumptions are made about the phase map in order to determine the dipole connections. However, if a series of two dimensional phase maps is available revealing the history of the appearance and annihilation of dipoles, the correct dipole connections can be directly observed without relying on any assumptions. In such a case, the third dimension will provide that extra information. 2.4.2 Three-Dimensional Poles To study three-dimensional (3D) poles, a 3D random phase map was produced. This phase map consists of a series of 2D phase maps varying gradually such that the first and last maps have no poles. The third dimension is considered to represent time.. For such a 3D map, 3D poles were detected and the presence of a pole (in any direction) was flagged. The result was a series of rings floating in three dimensions, clearly connecting the dipoles as shown in Fig.2.15. These poles, which are connected in three dimensions, would be separate in any 2D (x, y) slice. Clearly, future development of the pole-guided-cutline method of phase unwrapping presented in this chapter could consist of a 3D treatment of the poles. Such an approach is expected to overcome many of the limitations affecting 2D phase unwrapping methods. •3 x Figure 2.15: Resulting pole rings in three dimensions showing the dipole connections. 43 2.5 Summary In this chapter, phase maps have been characterised and certain structures within them defined. By understanding that the origin of the problem of phase unwrapping occurs due to poles, the source of these in MR images could be determined. This allows for certain assumptions to be made in order to best choose the path for unwrapping. Based on this, a new path-following, pole-guided-cutline method of phase unwrapping was developed and described in this chapter. Some results were shown as well as a test revealing that this technique is quite robust to the presence of noise. The novelty of the proposed phase unwrapping method lies in that it makes use of a new kind of quality map, called a score map. The score map is produced by repetitively tracing the dipoles connections, or fringelines, in the presence of several phase offsets. This score map, which indicates the relative phase gradient in regions of the image that are affected by poles, is quite unique. An advantage is that it is based only on the phase information, it does not rely on the magnitude. Furthermore, the development of the phase unwrapping technique proposed here is very general and hence there is hope of extending it to a three-dimensional method which could potentially overcome some of the limitations of two-dimensional phase unwrapping. 44 Chapter 3 The Use of Gradient Energy in Motional Ghosting Artifact Reduction 3.1 Effects of Mot ion in M R I The basics of 2DFT MR imaging was described in Chapter 1. Unfortunately, MRI is not as simple as it was presented there as is the case for most theories that apply to the physical world. The simplifications are crucial in allowing for the development of the theories but once the theories are to be put to practice, the complications of the real world set in. This thesis deals with the correction and better understanding of certain important artifacts that appear in MRI. In this chapter, the artifact dealt with is that due to motion of the object being imaged. In particular, it is the so-called ghosting artifact that will be the focus and this is a result of periodic or quasiperiodic motion. The theory for 2DFT MR imaging can be used to explain these artifacts that appear in real MR images of human bodies. Extensive investigations, both experimental and theoretical, have been done to understand the many effects motion will have on the resulting images [38, 39, 40]. Since these effects can often be detrimental, degrading image quality, masking or mimicking lesions etc., they affect the clinical interpretation and thus, diagnostic power, of MRI. It is for this reason that many researchers have attempted to correct or at least reduce them [41]. This is also the motivation behind the work presented in this chapter. A brief review of the theory explaining motion-induced artifacts in MRI along with an overview of some of the most prominent techniques used to deal with them are given in the first two sections of this chapter. The problem of motion-induced effects in MRI comes about because of a combination of two factors: a finite amount of time is required for MRI data acquisition and the human 45 body is not a still object. Although many fast imaging techniques have made MRI almost equivalent to a snap-shot in time [4], these techniques are still not widely used due to loss of signal-to-noise, spatial resolution and other image quality issues as well as hardware requirements. Still today, the standard sequences for imaging generally take about 2-30 minutes which means that any physiological motion occurring during this time is likely to have an effect on the resulting images. Although voluntary human motion (such as head or body shifts) can be somewhat controlled in most instances, the human body is a rich source of involuntary motion due to blood flow, cerebrospinal fluid pulsations, cardiac motion, gastrointestinal peristalsis and, to some degree, respiratory motion. Therefore, it is not surprising that some motion-induced effects are commonly observed in MRI. Furthermore, many of the involuntary motions are at least quasiperiodic in nature: breathing displaces the abdominal contents by 1 — 2 cm every 5 seconds on average while a cardiac cycle has a frequency of about 1 Hz [38]. It is the periodicity of these motions that leads to the ghosting artifact. Although ghosting artifacts depend on the exact motion that lead to them, the amplitude or intensity of the motion and the exact pulse sequence used to acquire the data, a general understanding of these artifacts is possible. The basic effects of quasiperiodic motion can be seen on an 2DFT MR image as blurring and the ghosting artifact which refers to spatially shifted replicas of the moving structures. Since the MRI signal is complex in nature, having real and imaginary parts, so too are the ghosts. This means that the ghosts may interfere with a background image by either increasing or reducing the overall magnitude and hence may appear as either bright or dark, depending on the other signal present in a given region of the image. It is clear that such an effect could greatly alter the appearance of the signal in some regions and hence could lead to misdiagnosis. Since TR, typically about Is, tends to be much longer than TE, typically about 50ms, it is likely that significant motion will occur during the time elapsed between the data acquisition and the next RF excitation pulse (approximately equal to TR — TE). The effects of motion occurring during this time interval will be referred to as interview effects as opposed to those occurring during TE which will be referred to as intraview effects. Most ghosting artifacts come from interview motion. In particular, breathing or other physiological motions which are relatively slow, will be effectively frozen in time with respect to TE and hence all of their effect comes about because of the motion during this TR — TE time. Faster motions (such as those associated with the cardiac cycle and blood flow) have an additional effect by causing a change in the imaged object during TE but this will not be further discussed here. In Chapter 4, which deals with the imaging of blood flow, the particular intraview effects of cardiac-related motion will be reviewed. Any interview motion can produce a change of signal intensity. Two such phenomena are 46 the "inflow" and "outflow" effects of flow into and out of a slice. The "inflow" effect is seen in gradient echo images where "fresh" incoming spins, that have not been excited before, can produce a stronger signal than the stationary spins which may not have completely recovered after TR. The "outflow" effect can be seen in spin echo images where the spins must experience both the initial RF excitation pulse and the 180° pulse in order to contribute to the signal, otherwise resulting in signal attenuation. These effects can be useful for flow visualization techniques [4]. Local blurring is a common detrimental effect of motion apparent in many imaging techniques where the signal produced is a time average. However, ghosting is a very particular kind of motion-related artifact, directly related to the use of Fourier transforms in MRI. Since the data acquired to produce an MR image is acquired in the Fourier transform space called &-space, the ghosting artifact can be explained by looking at how motion affects ft-space acquisition and how this effect is translated into the image space. As described in Chapter 1, an induction signal in MRI, kx, is acquired for one phase-encoding value, ky, at a time. After a time, TR, the procedure is repeated with a different value for ky and this repeats Ny times. An inverse Fourier transform of the complete (kx, ky) space then yields the image. For most purposes, MRI can be considered a linear system where the resulting image is the convolution of the imaged object, with a point-spread-function (PSF) [38]. Motion occurring between ky acquisitions, will result in a modulation in the ky direction which, when inverse Fourier transformed, will result in off-centered peaks for the PSF in the y direction. These off-centered peaks represent the harmonics of the motion meaning that they reflect the amount of signal, of a particular frequency, contributing to the overall motion. The peaks closest to the center indicate the first order (or fundamental) harmonic of the motion meaning that they reflect the smallest frequency present in the motion which is in fact equal to the frequency of the motion itself. The further away a peak is from the central one, the higher order the harmonic (or frequency) it represents which means that it represents the more rapidly varying components of the motion. Since the overall motion has a particular frequency, given by the fundamental harmonic, all other harmonics contained in the motion will be multiples of this one, hence the equal spacing of the peaks. The ghosts are basically the result of a convolution of this PSF with the average magnetization being imaged and hence occur along the phase-encoding direction. The central peak of the PSF represents the DC (direct current) component of the motion or the time-averaged modulation of the signal. For no motion, ideally, one would expect this to be a delta function but, usually, motion causes a broadening of this peak. In fact, the effect of motional averaging, seen as blurring in the image space, results in a broadening of the PSF peaks and this will be proportional to the amplitude of the motion as is observed in many other 47 imaging modalities [40]. The important observation which has been made experimentally and derived theoretically [38, 39], is that the effects of quasiperiodic motion will be the same in nature, regardless of the direction of the motion. Although the blurring will be in the direction of the motion, the ghosting will always appear as a shift in the phase encoding direction due to the actual motion-related modulation of the ky direction. This is true for all in-plane motion along (x, y), even though motion along the x direction causes a modulation of the signal frequencies while motion along the y direction causes a modulation of the signal phase. Although the pattern may be more complicated for perpendicular motion, along the z direction, if the centre of the oscillation does not correspond to the centre of the slice and depending on the shape of the slice selected, for practical purposes, the following conclusions can be made about motion-induced ghosts. Regardless of the direction of the motion, the overall intensity of the ghosting is dependent on the intensity of the moving structures while the relative intensity of the ghosts (representing the different harmonics of the motion) depends on the amplitude of the motion: the larger the motion displacements, the more intense the ghosts that are further away from the source of motion. However, the spacing between the ghosts, yo, does not depend at all on the amplitude of the motion but rather on its frequency, / , as given by: _ 2TT • / • TReff V G ~ ^AG-ty ( 3 ' 1 } where 7 is the gyromagnetic ratio for protons, TReff is the effective repetition time between RF excitation pulses (equal to TR • NSA if there is parallel averaging of the signal with fixed ky value, NSA times before going to the next ky value), A G is the phase-encoding gradient step size and ty is its duration. Since the field-of-view (FOV) in the y direction is given by: FOVy = 1.^..t (see Chapter 1, Eq.1.9), one can also write: V G f-TR.„ = *%U- (3.2) FOVy J ~~"KJJ T where T is the period of the motion. This simplification shows that the amount the ghosts are shifted along the phase-encoding direction relative to FOVy is equal to the proportion of the motional cycle occurring between subsequent ky values. Eq.3.1 also shows that the ghost spacing, yo, can be somewhat controlled by choosing appropriate values of TReff, AG and ty. This point will be further discussed in section 3.3.7. 3 .2 G h o s t R e d u c t i o n M e t h o d s : D e g h o s t i n g Ghost patterns consist of regularly spaced replicas of the moving structures and deghosting aims to reduce or eliminate the intensity of these ghosts from an image. Since the intensity 48 of the ghosts depends on the intensity of the moving structures, one way to reduce the effect of motional ghosting is to reduce the intensity of the moving structure in the image. This is an approach that can be used if the moving structure is not of interest. The intensity of an object can be attenuated by adjusting the contrast due to appropriate selection of TR and TE. Also, the choice of the imaging sequence itself can affect the signal intensity since, for example, blood flowing through a single slice will appear more intense in a gradient echo image than in a spin-echo image. Other intensity attenuating methods include spatial or spectral presaturation, short Ti inversion recovery (STIR) and surface-coil imaging [41]. As previously mentioned, another factor affecting the relative intensity of the ghosting is the amount of displacement or amplitude of the motion. This also affects the amount of blurring so controlling the amount of displacement would be desirable to reduce both motion-induced ghosting as well as the associated blurring. Restraints or voluntary breath-holding for instance can be ways of reducing these motional artifacts, however, these methods require patient cooperation and this is not always possible for younger patients or those in diseased states. Other ways of producing images where the amplitude of motion is effectively reduced is by gating, triggering, retrospective gating and post-processing by using a realistic model of the motion [41]. However, none of these methods has proven to be very reliable in all cases and hence, other options are still being investigated. The most successful and widely used approaches of motion artifact reduction seem to be those that attempt to reduce the effect of the periodicity of the motion without addressing the motional blurring effect. In other words, these methods are meant to cause a deghosting in the image. The following work described in this chapter falls into this category of methods. Other such methods include averaging, which increases data acquisition time, pseudogating, which requires a priori knowledge of the motion which must also be very periodic and ordered phase encoding (OPE), which is technically difficult to implement [41]. 3.3 Two-Point Deghosting for Superimposed Ghost Pat-terns 3.3.1 Introduction When two images are acquired with similar ghost patterns appearing in the same locations, these are referred to as superimposed ghost patterns. Since the term point is commonly used to refer to a phase encoded signal acquisition corresponding to a particular image, two-point deghosting relies on the input of two separate images to reduce or eliminate the intensity of the ghosts. Images with superimposed ghosts can be obtained by two-point interleaving 49 the exact same sequence. However, for such a case, the resulting ghosts will be almost identical (except for the small phase shift caused by TR, the sampling delay between the collection of the same line in k-space in the two images). A simple average will improve the signal-to-noise ratio but it will not significantly reduce the motion artifact. In fact, this is not an effective use of a second acquisition since little new information is gained. Such an averaging is equivalent to averaging several kx data acquisitions at a given fixed ky value before moving to the next ky of data collection and is referred to as parallel averaging. An alternative to this type of averaging is called serial averaging and has been suggested as being advantageous for ghost reduction. In serial averaging, a complete set of /c-space is collected several times in a row. Each set results in an image so there is no interleaving which means that the time between repeated data acquisition is increased. This increase in time is advantageous for ghost reduction since it will be more likely that the resulting images will have ghosts varying significantly in phase and hence more cancellation can be achieved by averaging. In general, if the number of signal averages is NSA, then a y/NSA reduction of ghosting intensities can be achieved. However, the increase in time is also a disadvantage to this method since it increases data acquisition time. The idea behind the Three-Point Ghost Phase Cancellation method [42] and later the Two-Point Interference method [43] was to make use of the phase of the MRI signal and acquire images that have superimposed ghost patterns that are phase rotated with respect to each other. Since the ghosting is a result of a periodic or quasi-periodic modulation occurring between successive readout lines, for a standard acquisition, the modulation will be along the phase encode direction. Assume this direction is the y direction, then the modulation affects A;-space along the ky direction producing a motion-related point-spread-function (PSF) with peaks along the y direction. In order to obtain the interleaved images with superimposed ghost patterns that are phase rotated with respect to each other, a relative phase ramp is introduced to the motion-related PSF that causes the ghosts. This is accomplished by introducing a shift of Aky between the acquisition of ky lines in the two interleaved images. This means that while one image is acquiring line ky in A>space, the other image is acquiring line ky + Aky. Therefore, one image acquisition proceeds from — ky-max to +ky_max while the other proceeds from —ky_max + Aky to +ky_max and then back again from —ky-max to — ky-max + Aky. The later A;-space acquisition is then reordered into the conventional format. Such a shift of /c-space acquisition followed by a reordering resulted in the desired phase ramp in the motion-related PSF, causing a phase rotation of the ghosts [42]. 50 3.3.2 Review of Existing Two-Point Interference Method For the Two-Point Interference method [43], superimposed ghosts were combined in such a way that the ghosts were made to destructively interfere, resulting in ghost cancellation. The formulation of this method was derived from the basic assumption that the ghosts were of equal magnitude and hence a simple rotation would be enough for complete cancellation. Therefore, this equal magnitude assumption resulted in the need to determine a single parameter, 9, representing the rotation. Starting from the decomposition of both interleaved images, I± and I2, into the time-averaged part, I0, and the ghosts, g\ and g 2 respectively, the exact form for the rotation angle can be derived. In the following, all symbols, representing vectorial quantities such as image I\ and ghost gi, can be considered complex numbers. Begin by writing the two images as: h = h + gi , , , , h = I0 + g2 = I0 + 9ieiB> 1 } where 0g is the ghost angle representing the rotation between two superimposed ghosts of equal order. The weighted sum, Iw, is given by: Iw = wi • h + w 2 • I2 (3.4) where w\ + W2 = 1 in order to preserve the identical time-averaged signal, I0. Writing Wi = w and w2 = 1 — w one gets: Iw = w-h + (l-w)- J 2 . (3.5) Jo can be solved from Eq.3.3 as a weighted sum: Jo = Iw = Ia + icot{0g/2) • Id (3.6) where Ia is the average image and Id is one half of the difference image: la = ( / l + / 2 ) / 2 , . h = ( / l - / 2 ) / 2 . { 6 J ) It would appear that once the ghost angle 9g is known at any location (for any given ghost order), the weighted sum that causes complete cancellation of the ghosts can easily be calculated according to 3.6. However, 9g is not known in general and since the images also contain noise, the problem of ghost cancellation is somewhat more complicated. Fur-thermore, 9g is different for different orders of ghosts so it must be determined regionally. In order to solve this problem, the gradient energy was introduced [43]. 51 The gradient energy (Eg) is defined in a region, R, of an image , I, with real and imag-inary parts given as Real[I\ and Imag[I] respectively, as the non-negative value resulting from the sum over all the pixels in R, of the magnitude squared of the gradient: In [43], Eg(IW) was used as a criterion to select the appropriate weight. The rationale for this is based on the fact that Eg measures the amount (quantity, size and magnitude) of edges in an image. The presence of motion artifacts is expected to increase the amount of edges in an image so when Eg is minimized, the amount of artifact should be minimized. By changing 9g to 9 in Eq.3.6, 9 was regionally varied until Eg(IW) was minimized in each region. This so-called "tuning" procedure resulted in good ghost cancellation and it was argued that in some regions where the noise was significant, such "tuning" would provide a good balance between ghost cancellation and noise amplification [43]. In fact, in noisy regions, the desired result is to perform a simple average hence Eq.3.6 shows that for this to happen, 9 = TT should be used, regardless of 9g. Eg minimization assures this. The result is that the single regional parameter, 9, required for best results is not always equal to 9g, it depends on the noise present in that region of the image. 3.3.3 Analytical One-Parameter Minimization Technique: l p G E M In [43], the "tuning" procedure involved stepping through a series of possible values for the parameter 9 and at each step, calculating the resulting IW as well as Eg(IW) and comparing the results in order to select the value of 9 that gives a minimal value of Eg(IW). This procedure was time consuming and it also made the quality of the result dependent on the step size. Due to time constraints, the value for 9 was determined regionally and then applied to all the pixels within the given region. This resulted in some sharp discontinuities between the regions in some areas of the result. Particularly, discontinuities were evident in the more homogeneous regions of IQ when both images were not identical there for some reason, see Fig.3.2. In this thesis, an analytical expression for Eg(IW) as a function of 9 is derived. From such an expression, an exact minimization of Eg(IW) with respect to the parameter 9 is possible. This results in a simple expression for the value of 6 that gives the wanted gradient energy minimization (GEM) and hence a single expression for the resulting Iw is obtained. Starting with the expression of IW as a function of 9 given by: IW = IA + icot(6/2) • ID, (3.9) 52 the real and imaginary parts of /„; can be written as: Real[Iw] = Real[Ia] - cot{9/2) • Imag[Id] . . Imag[Iw] = Imag[Ia] + cot{6/2) • Real[Id). [ ' Using Eq.3.8 and the above, an expression for Eg(Iw) as a function of 9 can easily be derived: Eg{Iw) = Eg(Ia) + cot2(9/2) • Eg(Id) + 2cot{9/2) • X{Ia, Id) (3.11) where X(Ia,Id) is a cross-term involving partial derivatives of Ia and Id and defined in a region, R, by: R dReal[Id] dlmag[la] dReal[Id] dlmag[la] dx dx dy dy dReal[Ia] dlmag[ld) dReal[Ia] dlmag[ld] dx dx dy dy (3.12) This analytical expression for Eg(Iw) as a function of 9 is a very important result. It not only insures that there is only a single minimum to the function that is being minimized, but now an analytical, exact minimization is possible. Setting the partial derivative of Eg(Iw) with respect to 9 equal to zero: dEg(Iw) d9 gives: „2 0, (3.13) csc^/2) [2cot(9/2) • Eg{Id) + 2X(Ia, Id)} = 0 (3.14) which results in: ° * W 2 ) = " W ( 3 1 5 ) for Eg(Id) ^ 0. This condition should be met everywhere except in those regions where both images are either identical or flat (contain no edges), although even there, the difference in the noise may be enough to assure this. Nonetheless, if this condition is not met, a simple average can be performed by setting 9 = TT whenever Eg(Id) = 0. The result for 9 given by Eq.3.15 can be used directly, with Eq.3.9, to calculate the value of Iw that has minimal Eg: ' - ' • - ' • ( w ) - ' * < 3 ' 1 6 ) The simplicity of this calculation allows for a sliding window application: X(Ia,Id) and Eg(Id) are calculated regionally but only applied to the central pixel. This is possible since the processing time has been dramatically improved by reducing the number of necessary calculations from 2TT/A9 to one, where A# is the "tuning" angle step size, in radians. Fur-thermore, no comparisons have to be made for this procedure. A typical value, yielding 53 reasonable results for the "tuning" procedure is A9 = 7r/30. This means that 30 calcula-tions and many comparisons have to be made for the "tuning" procedure whereas the new analytical implementation only requires a single calculation. 3.3.4 Comparative Results In order to assess the benefits of the improved, analytical implementation of the one-parameter GEM (lpGEM) technique for ghost suppression, it was applied to several images for which the original "tuning" was also applied. In all cases, the two applications resulted in qualitatively similar deghosting but the improvement in processing time was consider-able. An example of this is shown in Fig.3.1. where a stationary water bottle is imaged with a hand that is moving in quasiperiodic motion. The magnitude of one of the original ghosted images is shown along with the magnitude results of the two implementations of the lpGEM: a regional "tuning" and the analytical calculation applied in a sliding window. The only observable difference was that the images resulting from the sliding window ap-plication were smoother in regions where discontinuities had previously been seen. This is shown in more detail in Fig.3.2. where a close-up of (b) and (c) of Fig.3.1 are shown. (b) i **** (o) Figure 3.1: Comparison of magnitude images resulting from different implementations of the lpGEM. This example shows a stationary water bottle imaged with a hand that is moving in quasiperiodic motion, (a) shows the original ghosted image (Ji) (b) shows the regional "tuning" result for /,„ and (c) shows the analytical, sliding window result for Iw. 3.3.5 Discussion of Gradient Energy Curves According to Eq.3.11, the shape of the Eg(Iw) vs 9 curve is determined by two parts: the cot2(6/2) term which is symmetric and the cot(9/2) term which is asymmetric. If X(Ia, Id) is zero, only the symmetric part survives and Eg(Iw) will be a minimum for 9 — TT, resulting in 54 Figure 3.2: Close-up comparison of magnitude images resulting from different implementations of the l p G E M . (a) close-up of the regional "tuning" result for Iw showing the discontinuities (b) close-up of the analytical, sliding window result for /,„ showing no discontinuities. a simple average for /„, as seen by Eq.3.9. As long as X(Ia, Id) is not zero, another minimum will be found and the resulting weighted average will no longer be a simple average. It is expected that the cross-term X(Ia,Id) will be zero (or negligible) in ghost-free regions containing only the time-averaged signal. This is because each term in X(Ia,Id) involves a spatial partial derivative from both the average and difference images (see Eq.3.12). The difference image will be zero in these regions hence making one of the derivatives always zero so that each of the terms of X(Ia, Id) is canceled. In general, X(Ia, Id) is also negligible in regions of pure noise due to the high chance of having either Ia or Id being very small as well as the mutual cancellation of positive and negative terms in the summation. Since the behaviour of the lpGEM method depends on the location of the minima of the regional Eg(Iw) vs 9 curves, these curves were investigated. To study the behaviour of the regional Eg(I\y) vs 8 curves, they were plotted in several different regions of the image shown in Fig.3.1. As seen by Eq.3.11, these curves will be offset by a constant amount Eg(Ia) and this value may differ significantly in the different regions causing the curves to be spread over a large area of the y-axis. In order to make a better inter-region comparison, the curves were each normalized with respect to Eg(Ia) giving: This focuses on the behaviour of Eg(Iw) as a function of 9 only. Al l such normalized curves will pass through the value of Eg(Iw) = 1 when 9/2 = TV/2 and hence it will be easier to view them simultaneously. From Eq.3.11, it was observed that in the noise, where the two images are not correlated, the X(Ia, Id) term is negligible and hence the curves are symmetric (see Fig.3.3). In such regions, the lpGEM would result in performing a simple average since Eg(Iw) = l + cot2(8/2)-Eg(h) + 2cot(9/2) • X(IgJd) Eg(Ia) (3.17) 55 the minimum results at 9 = TT. In regions of ghosting, the curves are asymmetric, having a Sample F.g(Iw) v.v 0 curves for noise Trm 15000 0/2 • JT/2 -KXD0 5000 o ' ^^^^ 1 2 1 3 1.4: 1.5 i.e ti 1 -eooD -13000 0 /2 ( r a d i a n s ) Figure 3.3: i£<7 vs 0 curves in a region of noise. The curves are symmetric with minima occurring at 9/2 = TT/2. minimum on either side of 9/2 = TT/2 (see Fig.3.4). This minimum corresponds to 6g which is the angle between the two superimposed ghosts of same order. One would therefore expect the minimum to occur at 0™ = n • 01 where the n indicates the order of the ghost. However, upon further investigation, it was found that as the order of the ghost increases, the value of the minimum shifts away from the expected value, towards 9 = TT. This occurs since, as the ghost order increases, the signal decreases in intensity and hence the noise becomes more of a factor. The GEM criterion makes sure that noise amplification does not occur and a simple average or minimum at 9 = TT competes with the expected minimum at 9 = 9™. In the region where the discontinuities were observed, the curves were slightly asymmetric in all regions but largely varying minima occurred in each, thus causing the discontinuities. This was only observed in this example where I0 in the bottle was different in Ix and I2 although it was expected to be the same. A plot of several such curves taken in different regions of the bottle is given in Fig.3.5. An explanation for such a difference between the I0 in the bottle in the two images, has to do with the fact that the water present in the bottles may have undergone random motion which can result in significant phase differences in the interleaved images. If J 0 had been the same in both images as expected, Id would be zero and the calculation for 9 of Eq.3.15 could not be used. Instead, a simple average would be performed as stated earlier. This would not cause any discontinuities to appear. This is what happens in the I0 region of the hand. In summary, the Eg(Iw) vs 9 curves vary in symmetry according to the region. The result of this is that the minimum shifts along the 9 axis and thus the parameter determined by 56 Sample EgfJ^ vs 6 curves for regions of ghosts Ghost order +1 ,„ Ghost order +2 Ghost order +3 • -g£UUlJ •. ' , > . . 02 (radians) 0/2 (radians) 0/2 (radians) Figure 3.4: Eg vs 9 curves in a three regions of pure ghosts. Each set of curves correspond to regions within a given ghost order. It can be seen that within a single ghost order, if the signal is strong enough relative to the noise (as is the case for the first and second ghost orders) the minima of the curves fluctuate around a specific value for 9: 9/2 = 1.25, for the first ghost order, and 9/2 = 2.4 for the second ghost order. These values are close to having the expected relationship. The third ghost order region has lower signal and is more affected by noise as seen by the distribution of the curves' minima occurring around 9/2 = ir/2. the GEM criterion adapts in order to insure a balance between ghost suppression and noise amplification. In the noise, the curve is symmetric and the minimum corresponds to the parameter value which will lead to a simple average. In the ghosts, the curves are asymmetric (as long as the signal is strong relative to the noise) and the parameter corresponds to 8g so the result is a weighted average where the ghosts are made to destructively interfere. The curves are asymmetric whenever X(Ia,Id) is non-zero which can also occur in rare cases, when the IQ is not the same in I\ and I2- This can lead to discontinuities which are eliminated by the use of a sliding-window application of Eq.3.15. 3.3.6 New Generalized Two-Parameter Minimization Technique: 2pGEM In this section, a new generalized version of the above lpGEM method is presented. As stated in section 3.3.2, the single parameter that is to be determined for the lpGEM is a result of the equal magnitude assumption for the superimposed ghosts of same order. This led to the need to simply determine the rotation of the ghosts that would result in ghost cancellation. In order to generalize this method, the equal magnitude condition was relaxed and a two-parameter minimization, allowing for rotation and magnitude modification, was developed. The formulation of this method was developed starting from a new form for Eq.3.9. 57 DC W Sample F.g(Iw) v.v 6 curves for water bottle (I0) 50X00, 0/2 ( radians) Figure 3.5: Eg vs 0 curves in a region of Io- The curves' minima are distributed all over the place, no clear trend can be seen. This is what causes the regional discontinuities to appear when no sliding window is used. Decomposing Ia and Id in terms of the sought time averaged signal, I0, and the ghosts, gx and g2 gives: {a = I o + 9 a (3.18) h = 9d where ga = (gx + g2)/2 and gd — (gx — g2)/2. Eq.3.9 can then be written as: Iw = h + ga + i-w{9)-gd (3.19) where w(9) refers to a real function of the single parameter, 9. The exact form of w(9) is irrelevant. From this, another interpretation of the lpGEM is that a weight, w(9), must be found so that gd will cancel ga upon a rotation of 7r/2 (given by the i in front of w(9)), since this would result in Iw = I0. Therefore, the rotation is fixed at TT/2 and it is the magnitude of gd that will be modified by the parameter. Such a fixed rotation is a direct consequence of the equal amplitude assumption for the superimposed ghosts of same order. If gx and g2 are of equal magnitude, their difference will be perpendicular to their average. This is demonstrated in Fig.3.6 where gx and g2 are represented as vectors. Therefore, all that is needed for ghost cancellation is to rotate gd by 7r/2 and modify its magnitude so that it equals the magnitude of ga. This new interpretation led to a new, more general form of the GEM method for ghost suppression where the equal magnitude condition is no longer necessary. This form results from replacing w{6) in Eq.3.19 by a weight, w(A,B) dependent on two parameters, A and B, such that w(A, B) = A + iB: Iw = Io + ga + (A + iB) • gd. (3.20) 58 8il! - i;82 l a 1 ga Figure 3.6: Vector representation of gi and g2- In (a) the magnitudes are equal (||<7i|| = H52II) and hence the resulting ga and gd are perpendicular to one another. In (b) the magnitudes are not equal (||<7i|| ^  H52II) and hence the resulting ga and gd are not perpendicular to one another. These two parameters allow for any rotation of ga and magnitude modification in order to cause cancellation of ga. Using this expression and following the same procedure as before (Eqs. 3.8 and 3.10), one gets: Eg(Iw) = Eg(Ia) + (A2 + B2) • Eg(Id) + 2A • Y{Ia, Id) + 2B • X(Ia, Id) (3.21) where X(Ia, Id) is given in Eq.3.12 and Y(Ia, Id) is a similar cross-term involving the partial spatial derivatives of Ia and Id given by: From Eq.3.21, Eg{Iw) can be minimized with respect to A and B, by setting 3Eg(Iw)/dA = 0 and dEg(Iw)/dB = 0, yielding: Y(Ia,Id) = £ R + dReal[Ig] dReal[Id] dReal[Ia] dReal[Id\ dx dx dy dy dlmag[la] dlmag[ld] | dlmag[la) dlmag[ld] dx dx dy dy (3.22) 2A-Eg(Id) + 2-Y(Ia,Id) = 0 (3.23) and: 2B-Eg(Id) + 2-X(Ia,Id) = 0 (3.24) 59 which can be solved to give: A = -Y(Ia,Id)/Eg(Id) B = -X{IaJd)/Eg{Id). V-lb> The result of this two-parameter GEM method is that X(Ia,Id), Y(Ia,Id) and Eg(I,i) must be calculated in a region, R, and then the values can be used to determine the wanted Iw by substituting Eq.3.25 into Eq.3.20: -Y(ia,id) + i-x(ia,idy Iw — I a (3.26) Eg(Id) where Ia = Io + ga and Id = gd have been used. This will result in ghost cancellation for any superimposed ghosting pattern, regardless of whether the ghosts are of equal magnitude or not. This two-parameter GEM method (2pGEM) was tested on several of the images with equal magnitude superimposed ghosts and, as expected, the same results were obtained as for the lpGEM. The improved result for the 2pGEM was first seen when testing it on a simulated (synthetic) image with ghosts of differing magnitudes as shown in Fig.3.7. For this example, a central portion of the image is made to be equal in magnitude and phase in both images, I\ and I2. The upper and lower portions of I\ and I2 are made up of superimposed ghosts of differing phase and magnitude. Fig.3.8 demonstrates the improved deghosting resulting from the 2pGEM over the lpGEM on actual images with superimposed ghosts of differing magnitudes and phases. An explanation for how such images were obtained is given later, in section 3.6 3.3.7 Discussion on Window Size: Problem of Ghost Overlapping In the above, the value of Eg(Iw) has been minimized regionally meaning that Eg(Id), X(Ia,Id) and Y(Ia,Id) are all calculated for an n x n window and applied at the central pixel. However, no mention has been made about the determination of an appropriate window size yet the success of the 2pGEM relies on the use of an appropriate window. This is the case for any image processing technique that makes use of quantities that are inhomogeneously distributed over the entire image. In some cases, the result may be very dependent on the window size. In this section, the constraints on the window size for the 2pGEM method are discussed. If these constraints are satisfied, there is a range of values for the window size which will produce very similar results. The two parameters, A and B, required for the method of ghost suppression described above, depend on the value of the ghost rotation, 6g, the magnitude of the ghosts as well as the relative magnitude of the noise. Both 9g and the ghost magnitudes vary over the 60 (a) (b) (c) (d) Figure 3.7: Comparison of l p G E M and 2pGEM on simulated images, I\ and h, of superimposed ghost patterns with differing magnitude and phase. The magnitude images of I\ and I2 are shown in (a) and (b) respectively, with corresponding central y-profiles shown on the left and the phase of the ghosts shown on the right. The central portion of both images is identical, except for noise, and corresponds to IQ. The magnitude images of the results of the l p G E M and the 2pGEM are shown in (c) and (d) respectively, with corresponding y-profiles shown on the left. different ghost orders and hence image regions. There are no single values of A and B that will minimize Eg(IW) in every region. For this reason, the window cannot be the size of the entire image but rather, it should contain only one order of ghosting. The maximum size of the window is therefore constraint to be such that only one ghost order is covered at a time. This maximum window size constraint therefore depends on the spacing between the ghosts. This spacing is discussed in section 3.1 and Eq.3.1 shows its dependence on the frequency of the motion, / , and several parameters that can be somewhat controlled. The ghost spacing in number of pixels, PQ, can be written in terms of the ghost spacing, yc, given by Eq.3.1 and the field-of-view in the y direction, FOVy, as: P° = N"$yt (3-27) where N is the number of phase encode steps and hence number of pixels in the phase encode direction. Using this, Eq.3.2 and TRejf = 2 • TR for two interleaved images, one can write: P 0 = (2 • N • TR)/T. (3.28) In other words, the interghost spacing measured in number of pixels, PQ, is equal to the total number of cycles of motion sampled since the total scan time, T 0, is equal to 2 • N • TR. The interghost spacing, for a given FOVy, can be therefore be increased by increasing TR relative to T. There is, however, another constraint on this increase since as TR is increased relative to T, the quasiperiodic motion becomes less well sampled. If TR is too long compared to 61 *i ' '-4 ** - % (a) % V-6 I Cb) V (o) - * ' ' * V V % % * % • % % • i fife > * • v • *•% * r % • -, '••> % % Figure 3.8: Comparison of l p G E M and 2pGEM on actual images, Ix and I2, of superimposed ghost patterns with differing magnitude and phase. These images are the result of some processing as described in section 3.6. The experiment involved five tubes of water placed on a stand that is tilted back and forth in quasiperiodic motion. The magnitude images of I\ and J 2 are shown in (a) and (b) respectively. The magnitude images of the results of the l p G E M and the 2pGEM are shown in (c) and (d) respectively. T for the motion being imaged, the distinctly distributed ghosting pattern, representing the harmonics of the motion, tends to disappear and a more randomly distributed artifact results. This is somewhat related to the Nyquist sampling theorem which states that the highest possible frequency resolved is equal to half of the frequency of sampling. For periodic motion, any frequencies present in the motion higher than a half of the sampling frequency will appear as lower frequencies due to the phenomenon of "aliasing" (see Chapter 1, section 1.1.7). However, for the case of quasiperiodic motion, a randomness results if the motion is not sampled at a high enough frequency (short enough TR). This is often observed for cardiac-related motion, in images where the respiratory ghosts are nicely distributed, due to the much shorter period of the cardiac cycle relative to the respiratory cycle. These factors should be considered when selecting a particular value of TR for the scan. On the other side of the spectrum, the minimal value of the window size is usually con-strained by statistics. In other words, the larger the window size, the better the expected results due to averaging which tends to dampen the noise. However, since the GEM method results in a balance between ghost suppression and noise amplification, any local anomalies due to noise should be taken care of by the minimization and hence this may not be an im-portant consideration in this case. In fact, two other factors seem to play a more important role here. For the 2pGEM, it is the gradient energy that is determining the appropriate weighting factors for ghost suppression. This means that only the edges are contributing to the correct calculation while the flat regions (no spatial variation in magnitude or phase) favor simple averaging. If a window, centered on the background noise, is such that a bit of the ghost 62 edge is sampled, averaging will cause the ghost edge to contribute to the calculation, causing a shift away from simple averaging in that region. This means that averaging can cause "leakage" of the weights out to regions outside the ghosting. The same will occur as the weighting from outside the ghosts "leaks" into the ghosted regions. The result is a kind of blurred effect surrounding the ghosts and incomplete ghost cancellation within the ghosted regions as shown in the bottom row of Fig.3.9. This effect is small and is most noticeable for a large window size and in regions where the ghosts are surrounded by a flat background with no significant signal. Although, in general, this effect is not very significant, it does favor the use of the smallest possible window size for the regional calculations. In fact, the usual benefits of using a larger window and averaging do not seem to apply in this case since this "leakage" between different regions is undesirable. However, this is not always the case. In fact, such "leakage" is desirable when it occurs from the edges of the ghosted regions into a region within the ghosts that may be so flat that the GEM criterion does not work. If the ghosts of the images are very flat in a region and the window size is smaller than this region, no edges will be detected when the window is completely contained within the region and the 2pGEM will incorrectly reduce to simple averaging in these ghosted regions. The result is a minimal window size constraint: the window must be large enough such that, when centred on a ghosted region, an edge is always sampled. Fig.3.9 demonstrates the effect of different window sizes used to suppress flat ghosts for a synthetic image. If the window is too small, the central portion of the ghosts cannot be suppressed. If it is so large that it covers more than one ghost order, other portions of the ghosts cannot be suppressed. It is important to note the profiles on the left of all images. The effect of inter-region "leakage" is the slight blurred effect around and into the ghosted regions. It is noticeable in the images however, the profiles indicate the relative value of this effect is very small compared to the other effects such as incomplete cancellation of the central portion of the ghosts due to too small a window to sample the ghost edges (shown in (c)) and due to too large a window such ghost order overlap occurs (as shown by the arrows in (d) and (e)).. In summary, the window can be only as small as the image allows so that it always contains some ghost edges when centred on a ghosted region and only as large as the image allows so that only one ghost order is covered at a time. If the size of the flat regions in the image are not very significant, the window size can be made very small. This is the case for most actual images, although it does depend on resolution and field-of-view (FOV), it has been found that for most images, a window size of 3 x 3 is sufficiently large. Furthermore, for any window size of this value or greater, as long as it does not violate the maximum window constraint, the result is not very window dependent (see Figs.3.10 k, 3.11). 63 Figure 3.9: The effect of different window sizes used to suppress flat ghosts for synthetic images, (a) and (b) represent the magnitudes of the originally ghosted images with central portions corresponding to Io and the ghosts, above and below it, of differing phase and magnitude. Magnitude results of the 2pGEM are shown in (c)-(e). In (c), the window size is too small, in (e), the window size is too large and the leftover ghosting due to ghost order overlap is indicated by the arrows and in (d), the window size is both too small to suppress the central portion of the ghosts and too large so that ghost order overlapping is a problem. The window sizes used are indicated by the square box in images (c), (d) and (e) and vary from 19 x 19 to 35 x 35. All images have their central y-profiles shown along the left. Now that the constraints on window size have been discussed, the problem of ghost overlapping can be appreciated. In the case where different orders of ghosts overlap with each other, it is clear that there will be no window size small enough to satisfy the maximum window size constraint, i.e. any window will include more than one ghost order. However, it is not clear that this will necessarily cause an incorrect result. To further discuss this problem, the two images used for the deghosting will be considered to be made up of two superimposed parts: a background, made up of I0 (common to both images) and some noise (which is not common to both) and the ghosting patterns which are similarly distributed (same ghost overlapping occurs in both images) but which differ in magnitude and phase. What happens if two ghosts of different orders (rii and n2) completely overlap (edges are superimposed)? This could be the result of large interghost spacing and the wrap around 64 (a) it **K Figure 3.10: Effect of varying window sizes on actual images. These are images of a stationary bottle of water with a hand that is moving in quasiperiodic motion, (a) shows the magnitude image of one of the original, ghosted images, (b)-(e) show the magnitude results from using the 2pGEM method with increasingly larger window sizes: 3 x 3, 15 x 15, 45 x 45 and 75 x 75, as indicated by the outlined square boxes. effect where ghosts that should appear at line m appear at line m — N when m > N, where N is the image matrix size. In such a case, one can consider a combined ghost gni+n2 — dm + 9n2 which would be indistinguishable from gni. Therefore, the 2pGEM method should still be applicable in this long as there is still some difference (either in phase or in magnitude) between the corresponding gni+n2 m images I\ and I2. The problem of ghost overlapping arises when the ghosts do not overlap entirely. In these cases, there is a region that is made up of a combination of ghost orders that is surrounded by different single orders. Closer study of this situation revealed that the 2pGEM, which uses first-order spatial, partial derivatives, relies on the edges of the ghosts being surrounded by the background. The 2pGEM is successful in deghosting only when this is the case. For the fully overlapping ghosts where, in particular, the edges of ghosts overlap each other such that the resulting combined ghost still has all its edges surrounded by the background, the 2pGEM can be applied successfully. However, for the case of partial ghost overlap, the combined ghost no longer has all its edges surrounded by the background and hence, the 2pGEM is not as successful. This is demonstrated in Fig.3.12. 65 (a) (b) (c) (d) Figure 3.11: Effect of varying window sizes on actual images. In row I, the images are the cross-section of a chest and respiratory ghosting is present. In row II, the experiment involved 5 tubes of water on a stand that is tilted back and forth, in quasiperiodic motion, during the scan. Column (a) shows the magnitude image of one of the original, ghosted images, (b)-(d) show the magnitude results from using the 2pGEM method with increasingly larger window sizes: 3x3 , 15 x 15 and 35 x 35, as indicated by the outlined square boxes. Since the 2pGEM requires that the edges of the ghost be surrounded by the background for successful ghost suppression, some attempts of varying the window shape have been made. By varying the window shape, the hope is to find a window that only samples the edges of the combined ghost portion that are surrounded by the background. However, in general, it is not possible to find such a window for all regions of the combined ghost portion as shown in Fig.3.13 row II. There is, nonetheless, a particular case for which this solution of varying the window shape is very useful: when the ghosts are not actually overlapping but have a relatively large flat region relative to the ghost spacing such that they are effectively overlapping. When the minimal square window size must be so large that it samples more than one ghost order, the smallest and largest window constraints cannot be satisfied simultaneously and the ghosts can be considered as effectively overlapping. This is the case that was shown in Fig.3.9 (d), where the medium window size results in ghost order overlapping (as shown by the arrows) even though the window is still too small to successfully suppress the central portion of the ghost. Better ghost suppression can be obtained in these cases by varying the shape of the window so that the large dimension, 66 Figure 3.12: Effect of varying window size on synthetic images showing ghost overlap, (a) and (b) represent the magnitudes of two originally ghosted images with central portions corresponding to IQ and ghosts, above and below it, of differing phase and magnitude, (c)-(e) show the magnitude results from using 2pGEM with varying window sizes: 19 x 19, 29 x 29 and 35 x 35, as shown by the outlined square boxes on the right of each image. Al l images have their central y-profiles shown on the left. required to sample ghost edges when centred on the ghost, is perpendicular to the ghosting direction. An example of this is shown in Fig.3.13 row I. It is possible to imagine other ghost shapes and overlappings that would best be treated using more complex window shapes and sizes and maybe even a variation of these over the entire image. The result of the 2pGEM method can indeed be very sensitive to the window size when the constraints, described above, cannot be met. The unsuccessful results shown in Figs.3.12 and 3.13 suggest that maybe a combination of these results would give best results since, in many cases, the regions of unsuccessful ghost suppression seem to compliment each other. In other words, the regions with good ghost suppression for a small window size are, in general, different than the regions of good ghost suppression for a larger window size. However, simply re-running the 2pGEM method on two such complimentary results has not produced any improvement over simple averaging of these. This is probably because the results from the 2pGEM method have very smooth profiles as a result of Eg minimization and as such, the GEM criterion re-applied on these is not very effective. To 67 II Figure 3.13: Effect of varying the window shape on synthetic images. In row I, the ghosts are effectively overlapping if square windows are used. In row II, the ghosts are actually overlapping. Column (a) shows the magnitude of one of the originally ghosted images with central portion corresponding to Io and ghosts, above and below it. Column (b) shows the magnitude results from using 2pGEM with a square window of size 35 x 35 as outlined by the square box on the right of the images. Column (c) shows the magnitude results from using 2pGEM with a 1 x 39 window. Al l images have their central j/-profiles shown on the left. really benefit from combining the different 2pGEM results, another weighted sum could be used where, for example, a binary choice of one of the images is made at each pixel. For images such as the ones shown in these examples, an appropriate criterion for selecting the image at each pixel may be the absolute energy contained within a window. Another criterion may be more suitable for images with ghosting overlapping the time-averaged signal. Such considerations must still be investigated further but it seems likely that ghost suppression could be improved by finding appropriate combinations of window size, shape and results from the 2pGEM method. In conclusion, ghost overlapping presents a problem to the 2pGEM method and as of yet, there is no satisfactory solution. Fortunately, the ghost spacing can be somewhat controlled by the imaging sequence parameter TR and hence ghost overlapping can, in general, be avoided. 68 3.4 A New Application for 2 p G E M : Dissimilar Ghost Patterns 3.4.1 Introduction As stated previously, ghost patterns arise from quasiperiodic motion occurring between a given data acquisition line and the next RF excitation pulse causing inconsistencies between the several readout lines of A;-space. In general and the only case described thus far is when the readout direction is perpendicular to the phase encoding direction in /c-space. This results in ghosts along the phase encode direction because the motion modulates the signal between successive readout lines in £;-space and henceforth, in a direction perpendicular to that of the readout. This modulation leads to off-centered peaks of the PSF along the direction of modulation which results in the production of ghosts in the corresponding direction in image space. In general, the x direction is the readout direction and hence modulation affects things in a. direction perpendicular_this, along the..y direction which happens to be the phase encode direction. This perpendicular, independent relation between readout and phase encoding is not, however, a necessary condition for the production of 2DFT MR images. As long as the phase encoding and readout directions are not parallel, a given 2-dimensional region will be spatially encoded because A;-space will be sampled. Furthermore, the directions for the readout and phase encoding do not have to be vertical or horizontal but this is the most common and simple configuration mostly since the choice of x and y axis is arbitrary and hence is usually chosen to be along the readout and phase encode directions. In short, the ghosting resulting from a modulation of the signal in /e-space along a direction perpendicular to the readout direction, could theoretically appear in any direction, depending on the imaging sequence. In this section, a new oblique /c-space scan (Oi<"-scan [44]) is used to obtain images where the ghosting pattern no longer lies in a vertical or horizontal direction but rather along any oblique direction. It does this by acquiring &-space data along an oblique direction which is accomplished by an oblique readout gradient. An oblique readout gradient is the result of having some gradient component along both the x and y directions simultaneously. Minor modifications of a standard spin-warp sequence are necessary in order to obtain an OK-scan: a gradient along the y direction is added to the usual readout gradient along the x direction. This method is proposed in order to obtain two interleaved images with dissimilar ghosting patterns, ones that are not superimposed (occurring in the same locations), by varying the readout direction between image acquisitions, Fig.3.14. The combination of two interleaved images with dissimilar ghost patterns for artifact suppression has been previously proposed [45, 46]. Such dissimilar ghost patterns can sim-69 II ; 4\ v ^^^^^^^ Figure 3.14: Magnitude versions of two interleaved images with dissimilar ghosting patterns. Each column shows one of the two images. The first row shows cross-sectional images of the chest with respiratory ghosting in different oblique directions. The second row shows cross-sectional images involving two stationary bottles of water with a looped hose containing water that is being pumped periodically. The ghosts from the two portions of the hose loop (going into and out of the page) can be seen to go in different oblique directions in each image. ply be obtained by swapping the frequency and phase encoding directions between image acquisitions resulting in ghosting along the horizontal direction in one image and along the vertical direction in the other. Unfortunately, such a simple procedure will affect more than just the direction of the ghosting artifact. In fact, the off-resonance effects due to chemical shifts and any geometric distortion due to field inhomogeneities will also be different in the two images. This can lead to a misregistration problem and is therefore undesirable. Virtual frequency encoding [47] has been proposed as an alternative to simply swapping the frequency encoding and phase encoding directions. When the virtual frequency encoding direction is at 45°, the frequency and phase encoding axes can be swapped resulting in two images with perpendicular motional ghosting artifacts but no other unwanted, dissimilar dis-tortions. The OX-scan also produces two such images with identical off-resonance effects since the readout trajectories are Half-Cartesian (or Semi-Cartesian) [48]. This means that one of the directions (either x or y) is being acquired in the usual way, with each new data sampling point corresponding to the adjacent pixel, while the other is not. Off-resonance effects of such an acquisition scheme have been shown to result in a phase ramp along the 70 direction which is sampled with a cartesian trajectory [48]. Since this direction is the same for the two interleaved images acquired by an OK-scan, the off resonance effects in the two resulting images will be identical. Regardless of the method used to obtain the images, the following study applies to any two interleaved images with dissimilar ghost patterns. Since the Off-scan is easier to implement on a scanner with limited time flexibility and since it allows for a variety of ghosting directions, it is the method of choice for the following study. 3.4.2 Results The 2pGEM presented in section 3.3.6 causes ghost suppression by combining two images with superimposed ghost patterns that do not have to have equal magnitude. The formu-lation of this method did not lead to any restrictions on the exact magnitude of the ghosts and even magnitude zero should be possible. This means that if one image has a ghost in a certain location while the other image does not have any ghost component in that location, it should still be. able to cause ghost suppression. This allows for extension of-the 2pGEM. method to images with dissimilar ghost patterns. The 2pGEM method was therefore tested on several pairs of images obtained by an Off-scan, where the ghosting patterns are criss-crossed. This is obtained by switching the sign of the y-component of the readout gradient, Gy, while maintaining the same x-component, Gx, between image acquisitions. Fig.3.15 shows a sample of two such ft-space trajectories. In this example, the simultaneous readout gradients along the x and y directions have the relation: Gx = ±2 • Gy. In this case, the ghosting will appear at t an - 1 (.5)° on either side of the vertical or at a = 90° ± t a n " 1 (.5)° yielding a = 116.6° and a = 63.4° for images Ix and I2 respectively. Also, 7/8 th of the usual, rectangular, &-space will be sampled for the same number of phase encode steps. For two images, this represents 75% of the data that would be sampled the usual way. Fig.3.16 shows the effectiveness of the 2pGEM at suppressing dissimilar ghosting artifacts. 3.4.3 Magnitude Version of 2pGEM: General l p G E M This is the only section of the entire thesis which describes a method that does not rely on the phase information contained within an MR image. The method for ghost suppression described here results from a simplification of the application of the 2pGEM method previ-ously presented, for the special case of magnitude-only images. No novel ideas are presented in this section but some new insight into the behaviour of the 2pGEM method can be gained from the following formulation. Furthermore, by demonstrating that the 2pGEM method can be adapted to treat magnitude-only images, its applicability is expanded to many other 71 Figure 3.15: Representation of two interleaved, fc-space trajectories for the OK-sca.ii. In this example, the relation between x and y readout gradients: Gx and Gy, is: Gx = ±2 • Gy. The ghosting appears at tan-1 (.5)° on either side of the vertical and is defined in terms of the angles: ai = 90° + tan~1(.5)° and a2 — 90° - tan'1 (.5)°. In this case, 7/8th of the usual, rectangular data is sampled for each image, giving a total of 75% of the original data. areas of imaging. So far, the gradient energy has always been calculated in complex form which is given by Eq.3.8. However, the gradient energy of a magnitude image could be calculated as: (3.29) where the magnitude of image, I, is the usual: Mag[I] — (Real[I]2 + Imag[I]2)ll2. There-fore, given any complex image, I, there are two possible ways of determining the value of Eg (I): complex form or magnitude form. If the magnitude form is used, the entire for-mulation for the 2pGEM simplifies greatly resulting in an equivalent lpGEM which does not rely on any equal magnitude assumption, this will be referred to as a General lpGEM. This reduction from a 2-parameter minimization to a single parameter minimization can easily be explained. In the 2pGEM formulation, the two parameters are there to allow for a magnitude modification as well as a rotation of Id (or gd) such that it is made to cancel the ga component of Ia (see Eq.3.20). In the magnitude-only formulation, the rotation is no longer necessary due to the following. If a magnitude only version is used, Eq.3.20 becomes: Mag[Iw] =Ma + (A + i-B)-Md (3.30) 72 ( a ) Cb) ( c ) Figure 3.16: Magnitude images showing the effectiveness of 2pGEM on images with dissimilar ghosting patterns. This example shows a cross-sectional image of the chest with respiratory ghosting. The interleaved, ghosted images are shown in (a) and (b) and the result of 2pGEM is shown in (c). where Ma and Md represent the magnitude versions of the average image (Ia) and half the difference image (Id) given by: Ma = (Mx + M2)/2 Md = (Mi - M 2 ) /2 . V-6L) where Mx = Mag[Ii] and M2 = Mag[I2]. Given only the magnitudes of the quantities, one can choose to represent them as lying along the real axis with zero imaginary component. In other words, Ma and Md lie in the same direction and hence only a magnitude modification is required for cancellation upon summation. According to Eq.3.30, one would expect to get B = 0. This is in fact what happens since Eq.3.25 is a function of X(Ia,Id) which has the partial derivatives of an imaginary component multiplied in each term (see Eq.3.22). Since the imaginary components and derivatives are all zero, cancellation of each of the terms of X(Ia,Id) occurs, giving: B = 0. Substituting Eq.3.31 back into Eq.3.30 and using the result: B = 0, will give: „ r r , M i + M 2 Mx - M2 Mag[Iw] = 2 + A • g (3.32) which can be rewritten as: Mag[Iw] = wx- Mx+w2- M2 (3.33) with wx = 1 + f and w2 = 1 — ^ such that wx +w2 = 1. The 2pGEM applied to magnitude-only images is therefore equivalent to a standard, single parameter weighted sum problem with the added criterion of GEM used to select the weights. This method was applied to images shown in Fig.3.14 and the result is shown in Fig.3.17 (c). 73 3.4.4 New Use of Gradient Energy: Conditional 2pGEM The 2pGEM method for ghost suppression has been shown to be successful for ghost sup-pression given any two initial interleaved images which have either type of ghosting pattern (superimposed or dissimilar) and where the only difference between the images consists of their ghosting and noise. It can be seen in Eq.3.26 that if Id is zero, then Iw is a simple average. However, if Id ^ 0 and in particular, Eg(Id) ^ 0, then Iw will not be a simple average but rather some weighted sum that minimizes the gradient energy. Such a mini-mization will generally produce some signal cancellation in that region. This is desired in a region of ghosting but what if, for some reason, the time-averaged signal, 70, is slightly different in images I\ and I-p. In such a case, the 2pGEM would cause some unwanted signal cancellation in regions where Eg(Id) ^ 0. This discrepancy between the I0 in Ii and I2 will result, for example, if the interleaved images are flow sensitized and flow compensated (see Chapter 4). The second row of Fig.3.14 shows the magnitude images of two interleaved images that are flow sensitized and compen-sated. These images have non-identical time-averaged signals from the pump. This is not very evident in the magnitude images since the magnitude does not vary as much as the phase. If the previously described 2pGEM method is used in complex form on these images, some unwanted darkening occurs in the region of the time-averaged signal from the pump as shown in Fig.3.17 and indicated by the arrows. A magnitude version of the 2pGEM does improve the result significantly due to the smaller magnitude variation there but the result is a magnitude-only image. In some cases, the real and imaginary components may be re-quired and hence another approach must be found. Although no good solution has yet been found to solve this problem in all first step is to be able to detect the regions where this unwanted signal cancellation has occurred. The gradient energy has proven very useful in addressing this problem. The main purpose of describing the following test is to show yet another important use of the information provided by the Eg of an image. Hopefully, this can provide some motivation for further exploitation of this value in all areas of image processing. Within the range of statistical error, the gradient energy has been seen to obey the fol-lowing relationship for any two images, I\ and I2 that have uncorrelated edges or derivatives [43]: Eg[I1±I2) = Eg[I1] + Eg[I2] (3.34) For dissimilar ghosting patterns, the two interleaved images, I\ and I2, are each made up of a common time-averaged part, I0, and the respective ghosting components g\ and g2. Using 74 Fi gure 3.17: Magnitude results from applying the 2pGEM in complex form, in (b), and in magnitude form, in (c), to the two images shown in the second row of Fig.3.14. These interleaved images were acquired such that one is flow compensated while the other is flow sensitized (see Chapter 4). The magnitude of one of these images is shown in (a). The I0 signal from the pumped water in the looped hose is indicated by the arrows. Eq.3.7 and the fact that one can assume that ga and gd are uncorrelated to I0, we can write: Given Eq.3.34 and the fact that for dissimilar ghosting, gi can be assumed to be uncorrelated to 0 2 , one can see that Eg[ga] = Eg[gd] since ga = [gx + g2)/2 and gd = (gx - g2)/2. From this and the above one gets: Eg[Ia] — Eg[Id] = Eg[I0}. This result suggests that one could determine the expected minimum Eg[Iw], equal to Eg[Io] and estimated by Eg[Ia] — Eg[Id]. In other words, if one were to test the value of the resulting Eg[Iw], if Eg[Iw] < Eg[I0] then it can be assumed that too much cancellation has occurred. This test was used on the example of Fig.3.17. When it was found that Eg[Iw] < (Eg[Ia] — Eg[Id}), another approach was taken: a simple average of real and imaginary parts of I\ and I2 was done. This resulted in a complex result, with both real and imaginary parts, containing the ghost suppression provided by the 2pGEM method but where the unwanted darkening effect has been avoided. Such an implementation of the 2pGEM method is referred to as a conditional 2pGEM. The magnitude result of this implementation of a conditional 2pGEM is shown in Fig.3.18. The dark spots seen in (d) correspond to the regions of (b) where the 2pGEM method produced two much signal cancellation according to: Eg[Iw] < (Eg[Ia] — Eg[Id]) . These regions are no longer there in (e), where /„, is the result of the conditional 2pGEM method, shown in (c). This study shows that such an Eg test could be used to determine when the 2pGEM has produced a satisfactory result. When the result is not satisfactory, another option, such Eg[Ia] Eg[Id] Eg[h] + Eg[ga] Eg[gd}-(3.35) 75 Figure 3.18: Results of the Eg test. Images (a) and (b) are the same as those shown in Fig.3.17 where an unwanted darkening occurred in (b) relative to (a). Image (c) shows the magnitude of the complex result produced by using the conditional 2pGEM described above. The images shown in (d) and (e) are the results of performing : Eg[Iw] — (Eg[Ia] - Eg[Id}) where /„, is given by (b) and (c) respectively. The dark regions represent negative values. as simple averaging, could be used instead resulting in a conditional 2pGEM method. 3.5 Gradient Energy as a Measure of Image Quality 3.5.1 Introduction For all the previous discussions, the effectiveness of ghost suppression has been given qual-itatively. This is customary for ghost suppression algorithms since there is yet no standard measure of image quality which quantifies ghost suppression. Most new methods are simply qualitatively evaluated for their effectiveness. A quantitative measure of ghosting has great potential. It would result in a more systematic evaluation of ghost suppression techniques and could potentially automate and optimize these methods, improving their efficiencies and abilities. Some attempts have been made at quantifying the amount of ghosting in terms of ghost intensities in different regions [49] and signal-to-noise ratios [50]. However, it 76 is not clear that a decrease in signal intensity in a ghosted region is a good enough measure of ghost suppression ability. Any ghosting that overlays a given structure may produce a decrease in the signal intensity in which case, any reduction of the ghosting would result in an increase in the intensity. This reflects the difficulty with quantifying any amount of ghosting which is why ghost suppression algorithms are generally, qualitatively evaluated. An algorithm that corrects motion artifacts using a quantitative criterion has been pro-posed by others [51]. The entropy of an image is defined in their paper as a measure of the contrast of the image. Since motion often causes blurring and ghosting which spreads the signal intensity into regions that would otherwise be dark, it is argued that the smaller the entropy of an image, the higher the contrast, the more the dark pixels and hence the less the motion induced artifacts. In other words, the entropy of an image is suggested as a quantitative, relative measure of image quality indicating the amount of ghost suppression. Although this argument can be made for the ghosting that corrupts an otherwise dark region of the image, it does not apply again in the case of ghosts overlapping a given structure. In this case, more contrast or more dark pixels may not indicate less ghosting and hence the entropy cannot be expected to be a consistent measure of image quality in all regions. In this section, the gradient energy (Eg) is presented as a new, relative measure of image quality. The argument for it being a good measure, correlating with the amount of ghost suppression, has already been given in section 3.3.2. This argument is similar to that given for the entropy [51] but in this case, Eg is a measure of the edges, or differences in pixel intensities, rather actual pixel intensities. It is this basic difference that makes Eg a more global criterion that can be consistently argued for ghosting that occurs in any region of the image. Basically, ghosting is expected to increase the amount of pixel intensity differences (by either increasing or decreasing the signal intensity in some pixels) and hence increase Eg. Eg is hence proposed as a measure to quantify the relative amount of ghosting in an image: the greater the Eg, the more the amount of ghosting and conversely, the less the Eg, the more the amount of ghost suppression. This measure has already been extensively exploited by the ghost suppression methods developed in this thesis as well as others [42, 43, 50]. McGee et al. have recently shown that Eg correlates with observations as predicted above [52]. They performed a systematic study where several metrics (quantitative mea-sures) of an image were studied [52]. To do this, a comparison was made between qualitative observer ranking of several different images and the ranking of the same images, resulting directly from the value of the metric. The purpose was to find a quantitative measure, or metric, that best reflects qualitative evaluation of image quality as defined by the amount of motion induced artifact contained in the image. In this work, 24 metrics were applied to several images containing varying degrees of ghosting. Four observers scored the same 77 images and a comparison was made between the metric values and the observer scores in order to determine which metrics correlated the best with observer scores. The result was that the normalized gradient squared (which is a normalized version of Eg) and the entropy of the one-dimensional gradient along the phase-encoding direction were the two metrics that were significantly better correlated with observer scores than any of the others. Since Eg involves a much simpler calculation than the entropy metric (that uses.a logarithm) and the result was very similar for both these metrics, it does indeed seem to be an appropriate choice for quantitative evaluation of ghost suppression. In the case of the study done by McGee et al, the Eg was normalized. This is necessary in order to use the value of Eg to compare different original (ghosted) images or the resulting processed images of several different original images, since Eg is a relative measure. However, if Eg is used to com-pare the image quality of images resulting from different ghost suppression algorithms but applied to the same original images, the normalization is not required. In this section of the thesis, a study of Eg values of images is presented in order to demon-strate how well Eg correlates with the observable quality of the images. In doing so, a quantification of the qualitative descriptions given throughout the previous sections of this chapter is given as well as a quantitative evaluation of the effectiveness of the 2pGEM method relative to other methods of ghost suppression. In the first part of this study, the total and/or regional Eg is calculated for some of the images previously shown in this chapter. In particular, Eg values of the images of section 3.3 are used in order to quantify some of the statements made about the quality of the images. Since quantification allows for an unbiased evaluation of ghost suppression algorithms, the second part of this study consists of using Eg values to compare the results from the 2pGEM method with those of other ghost suppression techniques. In particular, the 2pGEM is used on some of the images from section 3.4 along with another two methods of ghost suppression: (i) a geometric average [46] and (ii) a weighted average [45]. For simplicity, magnitude-only images of I\ and J 2 are used and referred to as Mi and M 2 . For this reason, magnitude-only images result. The corresponding Eg values of the result of a simple average, Ma = (Mi + M 2 ) /2 are used as a reference. The geometric average consists of cross-multiplying the two originals such that the result, Mga, is given by: 3.5.2 Method (3.36) 78 This is expected to improve the result from that of performing a simple average, particularly in the regions where the ghosts contaminate the background noise. In these regions, the lack of correlation between the ghosting in M i and M 2 will produce more cancellation than if a simple average of M i and M 2 is taken. The weighted average, Mwa, as given by [45], has the usual form: Mwa = wl-Ml+w2-M2 (3.37) where w\ + w2 = 1. However, each weight is calculated as being inversely proportional to the relative amount of noise in the corresponding image (i.e. if M i has much more noise than M 2 in a particular region, w2 should be much greater than WI). Since noise will tend to increase the amount of edges present in an image and in doing so, increase Eg, the relative value of Eg can be used as an estimate of the relative amount of noise in each of the original images. In other words: Eg(M2) "'' Eg(Ml) + Eg(M2) , , Eg(Mi) W 2 ~ Eg (Mi) + Eg(M2)' Using Eqs.3.37 & 3.38, a particular version of the weighted average proposed in [45] can be applied. 3.5.3 Results & Discussion In section 3.3.3, an analytical form for the lpGEM method of ghost suppression was derived. The result is a more efficient implementation of the lpGEM method, allowing for a sliding window application. Fig.3.1 compares the result of such an application to the result of using the original step-wise "tuning" implementation. It is observed that both implementations produce qualitatively similar ghost suppression although the effect of both is definitely quite different in the stationary bottle. In order to isolate the ghost suppression ability of both implementations, a regional calculation of Eg was performed in each of the resulting images. The region used for the calculations is shown in Fig.3.19(a) and the resulting values are presented in Table 3.1. In this case, the quantitative measure given by the regional Eg is consistent with observation since the values resulting in each image are identical to within less than 1% which can be considered well within experimental error. In section 3.3.6, the lpGEM and 2pGEM methods are compared. Fig.3.7 compares the result of both methods used on two simulated images with ghosts of differing magnitudes. It is quite obvious that the 2pGEM method results in better ghost suppression. Eg cal-culations were performed for the total resulting images as well as for a region of ghosting 79 Section Figure used for Total Eg Regional Eg Figure showing Eg calculations region 3.3.3 3.1(b) 942 464 3.19(a) 3.1(c) 943 O i l 3.19(a) 3.3.6 3.7(c) 1 020 238 6 925 3.19(b) 3.7(d) 1 005 624 884 3.19(b) 3.3.6 3.8(c) 28 165 714 — 3.8(d) 32 164 718 — 3.3.7 3.9(cl)&(c2) 139 743 439 3.19(c) 3.9(dl)& (d2) 61 448 204 3.19(c) 3.9(el) & (e2) 63 976 612 3.19(c) 3.3.7 3.10(b) 757 166 3.19(a) 3.10(c) 960 286 3.19(a) 3.10(d) 1 031 104 3.19(a) 3.10(e) 1 582 226 3.19(a) 3.3.7 3.11 1(b) 200 626 672 3.11 1(c) 251 484 320 • • • 3.11 1(d) 253 616 592 * simple average (N/A) 255 055 328 3.3.7 3.11 11(b) 161 761 3.19(d) 3.11 11(c) 219 048 3.19(d) 3.11 11(d) 291 054 3.19(d) 3.3.7 3.13 1(b) 63 976 612 3.19(c) 3.13 1(c) 46 952 229 3.19(c) 3.3.7 3.13 11(b) 225 165 3.19(e) 3.13 11(c) 236 372 3.19(e) Table 3.1: Table indicating various Eg values. When a regional Eg calculation has been performed, the last column indicates the image showing the location of the region being used.*: the image is not shown anywhere, hence the N / A refers to "not applicable". within these, as shown in Fig.3.19(b). The calculations give relative values consistent with observation as presented in Table 3.1. The regional calculations reflect a much more signif-icant improvement of the 2pGEM relative to the lpGEM where Eg decreased by 87%. The regional Eg values are more significant when only a "small" region (in magnitude or area) of the image is expected to be affected by a ghost suppression algorithm. The averaging, caused by calculating Eg over the entire region, will dampen the effect of the improved ghost suppression, especially if the ghost magnitudes are quite insignificant relative to the time averaged signal. For the case of Fig.3.8, actual images are used to compare the result of the 2pGEM and lpGEM on ghosts of differing magnitudes. In this case, the total Eg is calculated for each of the resulting images since a significant portion of the original images 80 is ghosted. The effect of better ghost suppression seen to result from the 2pGEM is reflected in the values presented in Table 3.1 where Eg decreases by 12% from the lpGEM result to the 2pGEM result. Figure 3.19: Images indicating the regions used for the regional Eg calculations. In all cases, the magni-tude images of one of the originally, ghosted images is used to show the region. However, the actual images used for the calculations are indicated in Table 3.1. In section 3.3.7, the window size is discussed. Fig.3.9 demonstrates the constraints on the window size by showing the effect of different window sizes. In particular, too small a window size is used to show the remaining central portion of the ghosts as well as too large a window size which includes two ghost orders and causes an incomplete cancellation of the edges of the ghosts. A medium size window is also used to demonstrate a result that has both incomplete central ghost cancellation as well as incomplete ghost edge cancellation due to both effects. In order to quantitatively assess the differences in the resulting images due to the varying window sizes, regional Eg calculations were performed. The region used contains only ghosts in the original images as shown in Fig.3.19(c). The results, given in Table 3.1, reflect a significant decrease in Eg as the window size is varied from 19 x 19 to 29 x 29 and the amount of remaining central ghost portions is greatly reduced. However, 81 as the window size is further increased to 35 x 35, the Eg value actually increases. This shows the other effect, due to ghost overlap, taking over. In order to determine the window size that best balances both effects, the regional Eg was calculated for results occurring from window sizes of 31 x 31 and 33 x 33 as well. The quantitative values reveal that the best choice for window size is 31 x 31 if the minimal Eg can be used as a measure of best ghost suppression. This is consistent with qualitative evaluation, although the qualitative evaluation is more difficult to make since the images do not look very different and are hence not shown here. A qualitative evaluation would be very much less reliable in such a case, thus demonstrating the potential of a quantitative measure of the amount of ghost suppression. The effects of window size are also demonstrated on actual images as shown in Figs.3.10 & 3.11. In the case of Fig.3.10 a regional Eg calculation, in the region of ghosting as shown in Fig.3.19(a), was used to assess the ghost suppression differences of the various results due to different window size applications. The results are shown in Table 3.1 and are consistent with observation, confirming that the result is not as dependent on the window size for values from 3 x 3 to 45 x 45 where the increase in Eg is about 36%. As soon as the window size is very large relative to the ghost spacing, as it is for 75 x 75, the result is significantly worse since Eg more than doubles from 3 x 3 to 75 x 75. For case I of Fig.3.11, the ghosting due to respiration, is very closely spaced and is not as localized as it is for most of the other cases studies thus far. This case is interesting because there is a lot of ghost overlap occurring. However, the 2pGEM method is quite successful for a small window size as shown in the resulting images and reflected by the total Eg values presented in Table 3.1. A simple average of images h and I2 was also performed, in complex form, and the total Eg compared. The quantitative results suggest that as the window size gets large (in this case around 35 x 35) the result is no longer significantly different from a simple average. For each result of case II of Fig.3.11, the regional Eg is calculated in a region of ghosts as shown in Fig.3.19(d) and the results are presented in Table 3.1. In this case again, the Eg values correlate well with observation and confirm that result is not as dependent on the window size for values from 3 x 3 to 15 x 15 where the increase in Eg is about 35%. As soon as the window size is very large relative to the ghost spacing, as it is for 35 x 35, the result is significantly worse as it increases by 80% from 3 x 3 to 35 x 35. In section 3.3.7, the problem of ghost overlap is discussed. Fig.3.12 demonstrates the dif-ficulty presented by overlapping ghosts by showing some unsuccessful results of the 2pGEM. A variation of window shape is suggested as a solution for some cases such as that shown in Fig.3.9 where the ghosts are effectively overlapping. The regional Eg values calculated in 82 the ghosted region shown in Fig.3.19(c) and presented in Table 3.1, reflect the fact that a 1 x 39 window does improve the result for the case of effective ghost overlapping. However, for the other case of actual ghost overlapping, a regional Eg calculation in the region shown in Fig.3.19(e) suggests that the 1 x 39 window size gives a quantitatively worse result. These quantitative findings are consistent with the qualitative assessment as shown in Fig.3.13. In the second part of this study, Eg values are compared for the resulting images from different processing methods: • 2pGEM : MGEM • Geometric average : Mga • Weighted average : Mwa • Simple average : Ma. The results are shown in Figs.3.20 &: 3.21. Eg values were calculated both regionally and for the entire images and the values are presented in Table 3.2. Fig.3.20: (a) Mx (b) M2 (c) Ma (d) Mga (e) Mwa (f) MGEM Total Eg 10 352 795 10 210 166 8 465 217 8 444 064 7 853 902 7 532 168 Eg in ROh 2 344 3 062 1 323 1 158 886 893 Eg in ROI2 107 377 69 176 55 297 58 102 47 371 46 476 Fig.3.21: (a) Mx (b) M2 (c) Ma (d) Mga (e) Mwa (f) MGEM Total Eg 24 972 24 615 19 812 19 487 18 554 18 464 Eg in ROh 1 836 1 132 1 281 1 230 1 117 1 093 Eg in ROI2 423 869 321 201 128 120 Eg in ROh 829 682 389 370 315 314 Table 3.2: Table of Eg values for quantitative evaluation of ghost suppression methods. The total and several regional Eg values are shown down each column for the figures as indicated in the first and fifth row. The values are seen to correlate well with observation. The total Eg values followed the following relationship: Eg(Ma) > Eg(Mga) > Eg(Mwa) > Eg{MGEM) (3.39) 83 ROIi RO^ ROI, ROI2 ^ (a) \ J • »') ROI> (b) (c) RO^ R O l J • r / 4 \ ) ROI>—i • (e) R O I > - I ^ ( t) Figure 3.20: Magnitude images showing the results of applying several different processing techniques on the ghosted images of the first row of Fig.3.14. The ghosted, interleaved images are again shown in (a) and (b). The results of the several different processing techniques are shown in (c) M A , (d) MQEM, (e) M W A and (f) Mga- Two regions of interest (ROIs) for the regional Eg calculations of Table 3.2 are also outlined in each image. ROI\ samples a region that only contains noise and ghosting while ROI2 samples a region where the ghosting is overlapping the wanted signal I0. which can be used to argue that the 2pGEM results the most effective ghost suppression. As expected, M G A was similar to MA in tissue but showed better ghost suppression in the air. Although MWA was similar to M G E M in most cases, the 2pGEM seems to perform better for ghosts overlapping tissue edges as can best be seen by the quantitative results of Fig.3.21 in ROI\ presented in Table 3.2. In all cases, MQEM and MWA were better results than MGA and MA. 3.5.4 Conclusion To conclude this section, it can be stated that the Eg is a useful measure of the quality of the images. In particular, it is a relative measure that correlates well with the relative 84 R O I , ROL ROT, ROI, I ROIj * -ROI, ROI2 ROL it 1 ROI2 • ROT, • ROI. • ROI, (d) (e) (0 Figure 3.21: Magnitude images showing the results of applying several different processing techniques on the ghosted images of the second row of Fig.3.14. The ghosted, interleaved images are again shown in (a) and (b). The results of the several different processing techniques are shown in (c) M A , (d) MQEM, (e) M W A and (f) M G A . Three regions of interest (ROIs) for the regional Eg calculations of Table 3.2 are also outlined in each image. ROI\ samples a region containing an edge of the wanted signal, 7o and some ghosting. ROI2 samples a region that only contains noise and ghosting while ROI3 samples a region where the ghosting is overlapping the wanted signal, IQ, in a flat region of IQ. amount of ghosting seen in the images. The absolute value is not very significant since it depends on the size of the region for which Eg is calculated as well as the scaling of the particular image. However, in comparing several results from using either different methods (e.g. geometric averaging vs 2pGEM) or varying implementations of a single method (e.g. varying window size and shape) of ghost suppression on the same originally ghosted images, the values of Eg have been shown to be useful. These values have been used to quantify the observable differences and, in some cases, to distinguish between slight variations that are not easily detectable by pure observation. The hope is that such a measure be used to automate an unbiased selection of optimal ghost suppression. Furthermore, the Eg could be used as a measure of image quality extending into other areas of diagnostic imaging. 85 3.6 Ghost Suppression with One and a Half Excita-tions or Less 3.6.1 Motivation The GEM methods for ghost suppression presented thus far, are based on acquiring two full sets of data (or the equivalent of two excitations) in order to obtain a single ghost-free image. This led to the search for a more efficient way of acquiring enough data to obtain comparable results. In this section, a new idea is presented which would, in theory, reduce the number of required original, ghosted images, to the equivalent of one and a half excitations or less. A full description of the proposed method of ghost suppression using only one and a half excitations is presented below. It is then followed by a discussion on how to expand this idea to less than one and half excitations. Although some results show the effectiveness of the proposed efficient technique, it also has several limitations which are discussed below. . . 3.6.2 Theory In MRI, the usual way an image, 7, is obtained is by taking the Fourier transform of the k-space data Sk, which is the spatial-frequency representation of the image: I = FT{Sk). (3.40) If Sk is modulated by a function m such that: Sk.=m- Sk (3.41) then Fourier theory predicts that an equivalent operation in the Fourier transform space will be: /' = FT(Sk,) = FT(m) * FT(Sk) (3.42) where * indicates a convolution. This means that if the acquired data in fc-space is modu-lated by m, the resulting image, will be the image I convolved with FT(m) which can therefore be referred to as the point-spread-function (PSF) of the image acquisition: /' = PSF * / . (3.43) It is well appreciated that the acquisition of a full set of data for Sk is somewhat redun-dant [53]. Half-Fourier imaging (HFI) is one suggestion that has been proposed in order to significantly shorten the amount of time spent on image acquisition without sacrificing image resolution [53]. In this section, another approach is taken where every other line in 86 fc-space is sampled. This can be done by either sampling only the even /c-space lines or by sampling only the odd ones. These undersampling schemes will result in A;-space data that can be described by Eq.3.41 where m e v e n and m0dd are given by: m, 'even { 0 if odd /c-space line . . 1 if even fc-space line ^ ' ' or _ J 0 if even A;-space line , . m ° M ~ \ 1 if odd fc-space line. { i A b > For the following, the phase encode direction will be considered along the y-axis. The Fourier transforms of modulation functions given in Eqs.3.44&3.45 contain two delta func-tions: one at the center, reflecting the DC component of the modulation and another at the frequency of the modulation: y = N/2, where N is the number of phase encode steps. In particular, one can write: FT(modd) = PSF0dd = S y = 0 + 5y=N/2 • exp(m) ' ' ' where, PSFeven and PSFodd are related through multiplication of a phase ramp with gradient equal to TT/(N/2) since they are related through a spatial shift, of one pixel, in the reciprocal /c-space. Using Eq.3.46 in Eq.3.43, gives: I'even ~ PSFeven * I = (5y=Q + ^y=N/2) * I I'odd = PSF0dd * I = (Sy=o + by=N/2 • exp(m)) * /. (3.47) Using the fact that: and 6y=0 *I = I (3.48) Sy=N/2 * / = SHN/2[I] (3.49) where SHN/2[-] is a shift operator which takes the argument contained in the square brackets and shifts it by y = N/2, the predicted effect of such sampling schemes can be determined. The phase factor of exp(iTr), appearing in PSFodd, causes a rotation from SHN/2[I] to —SHN/2[I} and the final results can be written: I'even = I + SHN/2[I] I'odd = I-SHN/2[I] (3.50) In other words, an undersampling scheme consisting of taking only every other line in /c-space will result in an image that is made up of two superimposed components: the original image and itself shifted in space by y = N/2. Furthermore, if the sampling scheme 87 is such that only the odd lines are sampled, the shifted version of the image will also be phase rotated by TT (equal to multiplied by —1). The shifted versions of the original image will appear as ghosts (i.e. replicas) and will be referred to from now on as undersampling ghosts. Similar ghosts appear in a fast imaging technique called echo planar imaging (EPI) [54]. In this case, a phase error, resulting principally from time-reversal of the even and odd A;-space trajectories, results in a modulation of ft-space at the highest possible frequency, resulting in what are called N/2 ghosts. The reason this term will not be used here is because the undersampling ghosts will occur due to any undersampling scheme and will not always consist of a single one at N/2. The term, undersampling ghost is meant to reflect the more general case. An example demonstrating the ghosts resulting from undersampling the A;-space by taking only every other line, is shown in Fig.3.22. The observation that the (a) (b) (o) Figure 3.22: Figure demonstrating the undersampling ghosts that result from undersampling the fc-space data corresponding to image (a) and then Fourier transforming it back to image space, (b) is the resulting image if only the even /c-space lines are taken, (c) is the resulting image if only the odd fc-space lines are taken. Only the real part of each image is shown. undersampling ghosts resulting from Fourier transforming only the even lines of /c-space are out of phase with the undersampling ghosts resulting from Fourier transforming only the odd lines, led to the proposed efficient method of ghost suppression. 3.6.3 Method The proposed method consists of combining the 2pGEM method with the even and odd undersampling schemes mentioned above. Basically, three interleaved images are collected. Just as in section 3.3, each interleaved image contains the same time averaged signal, I0, 88 and some motional ghosts which are phase rotated with respect to each other. The motional ghosts are represented by gi, g2 and g3 as shown in Fig.3.23. The phase rotation is again Figure 3.23: Figure representing three interleaved images, I\, I2 and I3, broken up into their various components: a common part, I0, corresponding to the time-averaged signal, and three ghosts: gi, g2 and g3, respectively. These ghosts are phase rotated with respect to each other by arbitrary angles: a and /?, and they are not necessarily of equal magnitude. obtained by introducing a shift of Ak steps in the phase encoding order (see section 3.3). However, what makes this method efficient is that each one of the three interleaved images is undersampled, obtained by taking only every other line of fc-space. This means that each of the three images requires only N/2 phase encode steps (where the full set of data would require N steps) or the equivalent of .5 excitations, for a total of 1.5 excitations. The idea here is to use the fact that the undersampling ghosts, resulting from an even line acquisition, are out of phase with those, resulting from an odd line acquisition. More specifically, if two of the images, Ii-even and I3-even, result from even line sampling and one, h-odd, results from an odd line sampling scheme, they can be decomposed as follows: h-even = h + SHN/2[h] + 9i + SHN/2[gi] h-odd = h - SHN/2[I0] + g2 - SHN/2[g2] (3.51) h-even = 10 + SHN/2[I0] + # 3 + S-HN /2 [#3] where SHN/2 [•] again represents the shift operator. The odd line sampling scheme image can thus be used to cancel the undersampling ghosts of I0, written as SHNf2[I0] in Eq.3.51, in the images resulting from even A>space lines (since Io is identical in all three images). This results in two intermediate images, Iie2o and I2o3e, given by : he2o = h-even + h-odd = h + 9l+2 + S HN/2[g\-2] (3 52) hoZe — h —odd + h —even — -^ 0 + 02+3 + SHN/2 [#3_2] 89 where gi±j indicates gi ± gj. These two intermediate images have similar properties to the two initial images, I\ and I2, used in section 3.3 for the 2pGEM method. More specifically, they both contain a sought identical component, I0, and some differing motional ghosts. It is hence proposed that they be used by the 2pGEM method to extract the wanted, deghosted image, Jo-lt is important to note here that the ghosts in the intermediate images are not necessarily of equal magnitude since |<7i+2|| 1102+311 and H01-2I  7^  H03-2II, in general, even if ||<7,|| = H02I = 110311- This is why this method of efficient ghost suppression would not have been possible without the development of the 2pGEM where the equal magnitude condition on the motional ghosts has been relaxed. An example of one of these intermediate images is in fact shown in Fig.3.11 (II) since these images are used to test the 2pGEM method in section 3.3. Furthermore, more motional ghosts are present in these intermediate images than in the original, fully sampled ones, due to the added component of the undersampling ghosts of <7i_2 and g^-2 as shown in Fig.3.24. It is therefore expected that the problem of ghost overlap discussed in section 3.3.7 will be more of a problem for this efficient technique. (a) ** % rs * m *, (b) • • ; • % Figure 3.24: Magnitude images of an intermediate image, Iie3o, in (a) and the corresponding fully sampled image, I\, in (b). These images show the cross-section of a set of 5 tubes of water, placed on a stand that is tilted back and forth with quasiperiodic motion during the scan. The time-averaged signal lies along the diagonal indicated by the arrows. The other replicas of the tubes are the ghosts. 3.6.4 Results In practice, the implementation of the required imaging sequence that would produce h - e v e n , h-odd and h-even is not always easy. Fortunately, the proposed method could first be studied without having to acquire any new images, since the wanted A;-space lines (i.e. even or odd) 90 could be taken from already existing images. The results presented here are therefore part of a preliminary test, to see if the proposed method would work and to compare the results to those obtained from 2pGEM applied to the fully acquired data sets, h and I2. An actual implementation of the efficient method proposed has yet to be performed and the difficulties that must be overcome are discussed in the following section. A special case of the above method is to use Ii = I2 and hence a full /c-space of Ix is acquired and only the odd (or even) lines of another image, 7 3, are needed. This will result in a simplification of Eq.3.52 since Iie\0 = I\\ h = Io + gi ,„ ho3e = h+gi+3 +SHN/2[g3-i] These two images can also be used to do 2pGEM even though I\ has less ghosting than ho3e- Two images, I\ and 1$, have been acquired by interleaving the same sequence with a relative shift of Ak between their phase encoding steps. The result is two images with equal time-averaged signal, To, and phase rotated motional ghosts. These images were originally acquired to be used for the 2pGEM method proposed in section 3.3. For this study, only the odd lines of the corresponding fc-space of image 7 3 are used, giving the image I^-odd-Also, only the even lines of I\ are used to create the image Ii-even- An example of the real components of I\-even and h-0dd are shown in Fig.3.25 as well as their sum, 7 i e 3 o , in order to demonstrate the canceling out of the undersampling ghost of IQ. In Fig.3.26, the results of performing the 2pGEM on i i and I l e 3 o and on the two fully sampled original images, h and J 3 , are compared. Also, one of the original, fully sampled, ghosted images is shown as a reference. The results of this efficient ghost suppression 91 method, requiring only the equivalent of 1.5 excitations, can be seen to produce a significant amount of ghost suppression relative to the original, ghosted images. However, the result is not quite as good as if two full sets Ii and I3 are used as can be expected. Nonetheless, considering the time saving factor, the proposed technique may be an attractive option. (b) * * i t w Figure 3.26: Comparison of the results obtained by applying the 2pGEM method on two fully sampled images (requiring 2 excitations), in (b), and applying it on a single fully sampled image with another intermediate image (requiring 1.5 excitations), in (c). (a) shows one of the fully sampled, ghosted images. 3.6.5 Discussion A new efficient technique for ghost suppression has been proposed and some tests have been shown to produce satisfactory results. However, the actual implementation of this technique has yet to be performed since a rather complex imaging sequence would need to be developed. The sequence needs to acquire a full set of data for one image and only every other line of A;-space for the other interleaved image for which the phase encoding steps are collected with a shift of Ak. These considerations require extra effort for practical implementation on a clinical imager. However, the proposed technique is limited by more than just the practical aspects. The problem of ghost overlap is a considerable limitation here since so many ghosts are present in at least one of the intermediate images as shown in Fig.3.24. It is well appreciated that the 2pGEM method will not produce good results if the motional ghosts are overlapping (see section 3.3.7) and although this can be somewhat easily avoided for fully sampled Iy and 72, it becomes a much greater challenge for their undersampled counterparts. Fig.3.27 demonstrates the problem since two fully sampled images that have ghosts with enough 92 spacing to produce good ghost suppression when 2pGEM is used, cannot be well deghosted by 2pGEM when the equivalent 1.5 excitations are used, due to ghost overlapping. Figure 3.27: Figure demonstrating the problem of ghost overlapping for the efficient implementation of 2pGEM requiring only 1.5 excitations, (a) and (b) show the magnitudes of the fully sampled images, I\ and I2. (c) shows the magnitude of an intermediate image I\e20 where the ghosts can be considered effectively overlapping if a square processing window is used, (d) shows the result from 2pGEM applied to two fully sampled images I\ and I2 whereas (e) and (f) are the results of the 2pGEM performed on (a) and (c) with varying window shapes: 15 x 15 and 1 x 19 respectively. Regional Eg values are indicated on each image near the box outlining the region. The central y-profile of each image is shown along the left. Ignoring these limitations, it is theoretically possible to extend this method to produce even more efficient deghosting. In theory, the proposed technique could be extended to the equivalent of 1 + 1/n excitations for n > 2. This could be achieved in the following way. Every time the fc-space is sampled such that every nth line is taken, there are n possible ways of doing this. For example, there are two ways of sampling every other line, by taking either all the even lines or taking all the odd lines. One can write the n possible sampling 93 k\ + k\+n + ki+2n + &1+37Z • • • 62 + &2+n + k2+2n + &2+3n • • • 63 + h+n + &3+2n + &3+3n • • • (3.54) Sn = kn + + kn+2n + &rc+3ra • • • where £i indicates taking the ith line of &-space. Any time the A;-space of an image, I, is undersampled, the resulting image will be made up of I plus undersampling ghosts of I. The fewer the lines of /c-space taken (i.e. the larger the n), the more the number of undersampling ghosts of I resulting due to the more complex PSF associated with such a modulation. The images can therefore appear very corrupt if very few lines of /c-space are taken. However, summation of a complete set of undersampled images will eliminate the undersampling ghosts of 7 0 , as long as the undersampled images contain the same IQ. This is to be expected since summation of a complete set of undersampled versions of a single.image, J , should give.back the original, fully sampled image, I: £ > ? [ / ] = ! . (3-55) i=i where s™[7] indicates the image resulting by sampling I with the sampling scheme: s™. The idea for efficient ghost suppression is based on two different summations of complete sets of undersampled images so as to obtain two intermediate images, with no undersam-pling ghost of To but with several motional ghosts (both due to the original ones and their undersampling versions). Then, the 2pGEM is applied in order to extract the wanted I0. Since one of the complete sets can be a fully sampled image, only a single intermediate needs to be formed. In theory, this efficient method only requires a single, fully sampled image, Ii, along with a different image, I2, that can be undersampled using the ith sampling scheme from above: sf [I2\. The difference between i i and J 2 comes about by incorporating a shift of the phase encode steps so that the motional ghosts are phase rotated with respect to each other (see section 3.3). These acquisitions will be equivalent to 1 + 1/n excitations. An intermediate image, Iint, could be formed by adding the complete set of images: s™^] and Yjj^i s7j[h}- The 2pGEM could then be used to deghost Ix and Iint as long as there are differences in the ghosts appearing in Ii and those appearing in Iint. The fewer the amount of /c-space lines taken from I2 (the greater the n), the greater the required amount of dif-ference between the motional ghosts of I\ and I2: gi and g2. This is because the motional ghosts present in Jj n t , gint, that are used to cancel those in i i , are made up of a combination of <?i and g2 where the proportion of gi increases as n increases. From Eq.3.53, one can see that must be used to cancel g\. In this case, gint = 51+3 and since n = 2, gint is made schemes as: 94 up of an equal amount gi and g3. However, in general, gint can be written as: 9 m = { ^ 1 ) ' 9 i + 0 ' 9 2 - ( 3 - 5 6 ) If gint is going to be used to cancel g\ by the 2pGEM method, they must be significantly different in phase and/or magnitude. Using the above expression for gint, one gets: 9\ - g i n t = ~'(9i~ 92) (3.57) which shows that the necessary difference between gi and g2 required for the success of the efficient 2pGEM implementation requiring only 1 + 1/n excitations, is about a factor of n more than that required for the regular implementation using 2 full excitations. This idea was tested on some simulated images shown in Fig.3.28. The method was performed for 1.25 excitations where a fully sampled image I\ was combined with an un-dersampled image where every 4th line in fc-space, of I2, was taken, i i and I2 represent interleaved images acquired with a phase encode shift since both have equal I0 but their ghosts differ significantly in both phase and magnitude. In order to avoid ghost overlapping, only a single ghost was used since sampling every 4th line in A>space produces three, equally spaced, undersampling ghosts. The results, as shown in Fig.3.28, demonstrate that good ghost suppression was possible with only 1.25 excitations. In fact, the result was almost indistinguishable from the result of using both fully sampled I\ and I2. Quantitative eval-uation revealed a slightly smaller regional Eg value for the result using 2 full excitations. However, the discrepancy is of same magnitude as that resulting from a regional Eg calcu-lation in an equally sized ROI within the noise only. This indicates that the better value for Eg is probably a result of improved signal-to-noise due to the full 2 excitations compared to the 1.25 excitations rather than an improved ghost suppression. Although this test yielded very positive results, it was not possible to find an actual image for which a similar procedure could be tested. The limitations of ghost spacing caused by the problem of ghost overlapping and the required amount of difference between the motional ghosts in I\ and J 2 makes this method, as of yet, not practical. There is some hope that the limitation, due to the necessary amount of difference between the motional ghosts in I\ and I2, may be overcome. The intermediate image produced in the above method uses an equal weighting of all the fc-space lines. If the sampling scheme for I2 involves only every nth line of /c-space, then it is combined with the complimentary sampling of A;-space lines for Ix such that the motional ghosts from I2 will contribute (l/n)th to the resulting ghosts in the intermediate image. However, there is the possibility that variable weights could be used for the combination of these complimentary /c-space lines. 95 Figure 3.28: Results demonstrating the implementation of the 2pGEM on the equivalent of only 1.25 excitations, (a) and (b) show the simulated images representing the fully sampled I\ and 72 respectively, (c) shows the intermediate image produced by adding Yllzt st [h] a n d sf[/ 2]. The results of the 2pGEM applied to two fully acquired images and the equivalent of only 1.25 excitations are compared in (d) and (e): in (d), the 2pGEM is used on (a) and (b) with a square window of 15 x 15 and, in (e), the 2pGEM is used on (a) and (c). A regional Eg value is indicated in each image near the box outlining the region. The central y-profile of each image is shown along the left. The undersampled images written: Ylj^i^h and s™I2, could be combined by weights w\ and w2 to give the intermediate image Iint such that: lint = (wi • £ s]h j + (w2 • S?I2) (3.58) with w\ + w2 = 1 so that the time-averaged signal, equal in I\ and I2, is not affected. This may allow for a greater variation in the resulting intermediate ghosts: g2„( and the undersampling ghosts of these. This would also affect the signal-to-noise ratio (SNR) but maybe a compromise between the two is possible. Other than the variation of the weighting of each of the undersampled images, there is another degree of freedom associated with the creation of the intermediate image. It lies in 96 the nature of the sampling scheme used for the second acquisition. What if the &-space lines of I2 were to be taken at random, in a scrambled fashion rather than the regularly "every nth line" sampling scheme described above? The resulting ghosts would again be affected and as long as the missing fc-space lines from I2 are taken from the fully sampled I\, the undersampling ghost of 7 0 can be eliminated from Iint. These ideas could allow for amplification of any small difference between the motional ghosts of Ii and I2 although their benefits may be greatly limited by SNR considerations. These possibilities do however merit further investigation and may lead to ways of over-coming at least one of the problems associated with the efficient ghost suppression method described in this section. However, until a satisfactory solution to the ghost overlapping problem is found, the use of the efficient ghost suppression described in this section remains greatly limited. 3 . 7 Summary This chapter has presented several different techniques dealing with the elimination or at least reduction of the motional ghosting artifact. The techniques presented are all guided by the criterion of gradient energy minimization (GEM) for which an analytical expression has been derived. The two-parameter GEM method of ghost suppression (2pGEM) has been shown to be of the most general form. The only requirement for its successful application is that the two interleaved input images have ghosts that vary significantly. This method has been successfully applied to the case of superimposed ghosts as well as for dissimilar ghost patterns. The specific case of the application of 2pGEM to magnitude-only images, resulted in an equivalent general single parameter GEM (lpGEM) method. This is the only section of the entire thesis that does not depend directly on the phase of the signal. Nonetheless, its derivation stems from a simplification of the method of 2pGEM which was developed by considering the complex nature of the signal. The conditional 2pGEM method was also presented as a special case of the 2pGEM method. It addresses the case when the assumption made by the 2pGEM method, of equal time-averaged signal in both input images, does not apply. Throughout this chapter, the gradient energy of the images has been extensively ex-ploited. This allows for a better understanding of possible uses for such a measure as was shown for the case for the conditional 2pGEM method where gradient energy relations were used to indicate problematic areas in the results. Finally, the ideas for the development of a method that uses less than one and a half 97 excitations (full data acquisitions) was presented and some simulations demonstrate its possible ability. All the methods discussed in this chapter are based on the application through a sliding-window which means that they all are subject to the limitations of window size and shape dependence as well as the problem of ghost overlapping. 98 Chapter 4 Fast Dynamic Flow Imaging: Making Use of the Motional Ghosting Artifacts and Phase Information 4.1 Introduction In the previous two chapters, specific artifacts occurring in MRI were studied and the focus was then to eliminate, correct or suppress them. In this chapter, a different focus is taken in that the understanding of an artifact is exploited in order to extract useful information contained in it. More specifically, a procedure for fast dynamic flow imaging is presented. Dynamic flow imaging refers to the process of extracting information about the time evolution of the average velocity of the spins contained within a certain region of an image. The ideas and theoretical developments presented in Chapters 2 and 3 will be used here since the idea for fast dynamic flow imaging makes use of velocity encoded phase maps as well as the motional ghosting artifact. The ghosting artifact, described in the previous chapter, is a result of interview motion (occurring between data acquisition and the next RF excitation pulse). The ghosting arti-fact is characterized by misplaced signal intensity from one region to other equally spaced regions along a direction that is perpendicular to the readout direction. The goal of ghost suppression techniques is to reduce the intensity of the signal contained in the ghosts. In this section, the intraview effect of motion (occurring during TE) is presented. Al -though any periodicity in the motion occurring during this time interval will also contribute to the ghosting artifact, a more important characteristic of the associated artifact is that it affects the phase of the acquired signal. In contrast to the ghosting artifact, this effect on the phase can be desirable since it allows the MRI signal to encode flow information. First, a brief review of the theory behind velocity phase encoding in MRI is presented. 99 This is followed by a review of some phase methods of flow imaging, specifically, those for velocity phase mapping. Using the understanding of these velocity phase maps combined with an understanding of the ghosting artifact presented in the previous chapter, a new method of dynamic flow imaging is presented. Dynamic flow imaging refers to the production of a time-dependent movie of the flow. Such a movie, which can be clinically useful, is presently available by so-called cine MR imaging [37]. The technique for cine imaging involves retrospective sorting and interpolating of the data in order to obtain full fc-space data sets for various motion phases which are then further Fourier transformed in order to obtain the time evolution of the signal. Such a strategy combined with cardiac gated or triggered flow imaging can be used to obtain movies of flow [37]. However, cine imaging requires reliable motion monitoring or triggering and results in long scan times. The hope is that the method for dynamic flow imaging, proposed in this thesis, will provide the same information that is now available from cine flow imaging but in a much shorter amount of time and without the high demands on hardware and patient preparation. 4.2 Theory 4.2.1 Review of Velocity Encoded M R I Phase Maps Intraview motion, occurring during TE, induces phase shifts [41]. TE is usually relatively short compared to some physiological motions, such as breathing, for which this phase shift effect will be negligible. However, for more rapid motions, such as those-related to the cardiac cycle and blood flow, this effect will be significant since the amount of motion during TE will be significant. In the following, the motion of importance can be assumed to be the flow of blood or cerebrospinal fluid (CSF). Other motions of fluids that affect the MRI signal are diffusion, caused by random Brownian motion of molecules and characterised by its incoherence, and perfusion, which is motion on the capillary level that is characterised by its coherence [55]. However, it is the macroscopic flow effects that are the focus here. In 1982, Moran was the first to suggest the use of bipolar gradient pulses in order to introduce velocity phase-encoding to the data [56]. His work was the first, general, theoretical description of quantitative velocity mapping. In 1984, velocity phase-mapping methods were successfully demonstrated [57, 58] which extract the velocity of the signal from the phase information. The basic ideas for velocity phase-mapping in MRI are presented here. As previously stated, if there is motion during the application of spatial encoding gradi-ents, the phase of the MRI signal can be affected. More specifically, the relationship between the phase of the signal, <f>, and the time-dependent, three-dimensional position of the spins, 100 « ( t ) A G • G Figure 4.1: Bipolar gradient waveform used to velocity encode MR images. In particular, the first moment is given by: Gto{to + A) and is independent of the time origin of the gradient. r(t) is given by: rTE <j>{TE) = 7 / r(t) • G(t)dt (4.1) Jo where 7 is the gyromagnetic ratio and G(t) is the time-dependent gradient waveform in all three directions'of space. When motion is present, r(t) is not a constant in time and it can be expanded in the following Taylor series: r(t) = T0 + v 0 t + ^a0t2 + ... (4.2) where r 0 , v 0 and a 0 represent time-independent position, velocity and acceleration vectors respectively. Since each component of r(t) -G(t) contributes separately to the phase cp(TE), consider the component r(t)G(t) along a single direction. 4>(TE) can be written in terms of this component, using the expansion of Eq.4.2, as: / RTE RTE 1 FTE \ (f>(TE) = >y lr0 j G(t)dt + v0JQ tG(t)dt+-aoJ t2G(t)dt + ...) (4.3) which are the zeroth, first, second, etc. moments of G(t). If the spins being imaged have a constant velocity, VQ, then if the first moment of G(t) is non-zero, Eq.4.3 implies that there will be a non-zero contribution to the phase of the MRI signal. This is the basis of velocity phase mapping techniques. In particular, a bipolar gradient waveform is used as shown in Fig.4.1 since its zeroth moment vanishes while its first moment is non-zero. For most flow, the velocity of the spins can be assumed to be constant during the relatively short time TE 101 and hence only the first two terms of Eq.4.3 need be considered. Motion due to flow can therefore be sensitized along any direction by applying a bipolar field gradient along that direction. Phase-contrast angiography makes use of this phase sensitization and compares the phase in two images acquired with and without the insertion of a bipolar gradient. Any phase differences are indicative of motion of the spins and is used to extract the distribution of moving spins which is then displayed as a map [4]. Velocity sensitized and compensated phase encoded images have been used to validate the technique of velocity mapping in vivo [59]. The basis for the velocity sensitization or compensation ability is based on the first two terms of Eq.4.3. If the first moment of G(t) is zero, then 4>(TE) will be the same for spins moving and those that are stationary yielding a velocity compensated image. On the other hand, if the first moment of G(t) is non-zero, the spins moving will gain an extra mount of phase and the <f>(TE) term will differ between these, yielding a velocity sensitized image. The phase difference between these two images, provided it is not wrapped (see Chapter 2), will then be a good representation of the velocity field map since it will be proportional to it. This phase difference operation is necessary since the absolute phase of the signal has other factors than (f>(TE), given by Eq.4.3, contributing to it. A phase difference operation can be used to eliminate these other factors, isolating the effect of velocity-induced phase shift only. In practice, the phase of the MRI signal will also be affected by factors such as the chemical shift, acquisition window off-centering and field inhomogeneities. Furthermore, the ambiguity in the position for zero phase can be disregarded if a phase difference is used instead. In the description of 2DFT MR imaging in Chapter 1, slice selection and frequency encoding were obtained by using bipolar field gradients (for spin refocusing purposes) along the 2 and x directions respectively. This means that the sequence shown in Chapter 1 is inherently sensitized to flow occurring along either the x or z direction. In practice, the waveforms used for slice selection and frequency encoding are made to be compensated for flow sensitization by a technique called gradient moment nulling (GMN) [41]. This technique usually involves the application of another pulse to the gradient waveform so that the first moment is made to vanish at the time TE. The main drawback of this technique is that it prolongs the minimum TE that can be achieved. Fig.4.2 shows a sketch of possible gradient waveforms for pulse sequences leading to velocity sensitized and velocity compensated phase images. For simplicity, the sensitization is introduced in the z direction and it is assumed that there is no flow in the x and y directions since these directions do not have flow compensated gradient waveforms. In velocity phase mapping, the phase difference between velocity sensitized and velocity compensated images is displayed. As previously alluded to, an important consideration, 102 c o m p RF a. f l o w compensation n ^ selection 1 flowsensitization RF OAT-scan readout period readout period Figure 4.2: Sketch of possible gradient waveforms for pulse sequences leading to velocity sensitized and velocity compensated phase images. For simplicity, the sensitization is introduced in the z direction and hence in the slice selection pulse.' It is assumed that there is no flow in the x and y directions since these directions do not have flow compensated gradient waveforms. 103 affecting the success of this method, is the fact that it is difficult to unambiguously assign the correct value to phase differences lying outside the range [—TT,7T) (see Chapter 2). For this reason, the sensitivity of velocity phase encoding is often compromised. The method of phase unwrapping presented in Chapter 2 can hopefully be used to overcome this limitation for the dynamic flow imaging proposed in this chapter. 4.2.2 Review of Dynamic Image Reconstruction: DIR The dynamic image reconstruction (DIR), proposed by Xiang et al. [60], demonstrated that the motional ghost artifact, described in Chapter 3, carries useful dynamic information and should not be discarded. As already described, the ghosts result from a convolution of the time averaged signal with the PSF. Off-centred peaks in the PSF arise when there is a quasiperiodic modulation of the signal. These peaks are a direct result of the periodicity of the modulation and represent the contribution, of phase and magnitude, from each of the harmonics of the motion. In other words, the information contained in the ghosts reflects the information contained in the PSF which describes the motion. In DIR, a summation of the ghosts, weighted by appropriate temporal phase factors, yields the time-dependence of the structures imaged and is therefore used to produce a time-dependent movie of the image. The basic equation used for DIR is as follows: p(x, y; t) = I0{x, y) + g(x, y + dn) • exp(inuj0t) (4.4) n where p(x,y;t) is the reconstructed time-varying image, I0(x,y) is the ghost-free image representing the DC component of the spin density, g(x,y) is the ghost mask and t is the time in which the spin-density distribution changes. The summation is over all ghost orders present in the image (n goes from — n m a x to + n m a x ) . The incremental step, dn, is used to obtain the contribution from each of the ghost orders. It is given by dn = nP0, where P 0 is the spacing between ghosts, in pixels, along the phase encoding direction and given by Eq.3.28. For this equation, the phase encoding direction is along the y direction and the ghosting is also assumed to occur along this direction since the readout was along the x direction for DIR. The time-dependence is given by the phase factors: exp(inu>ot) which describe the time evolution of the harmonics with LOQ being the fundamental frequency of the motion (not to be confused with the Larmor precession frequency). The value of UJQ does not have to be known a priori since it can be calculated from P 0 . Eq.3.28 gives the relation between P 0 , the total scan time, T 0 = 2NTR, and the period of the motion, T. Since ui0 = 2TT/T, it can 104 be calculated given P 0 and To: wo = . (4.5) -to Therefore, in order to apply this DIR technique for a given pixel, (x0, yo), within the region of the time-averaged signal which is the source of the motion, the value of P0 must be determined and the real and imaginary components of Io{x0, Vo) and g(x0, yo + dn) must be known. In general, Io{xo,Vo) can be taken directly from either of the interleaved images if the region is specified. However, some of the ghosts usually overlap the stationary signal, hence, the ghost mask must somehow be separated from the stationary signal in order to obtain g(x0, yo + dn). In DIR [60], three interleaved images with superimposed ghosts were used. Special processing of the data was then developed in order to extract the wanted ghost mask. It is important to note that, in this case, PQ was used to step between ghost orders as well as to calculate the fundamental frequency of the motion, cu0- This is only a consequence of the fact that the ghosting direction coincides with the phase encode direction for DIR. This will not be the case for the oblique ghosting resulting from the OK-scan, as will be described for the proposed method. 4.3 Method 4.3.1 Proposed Method for D y n a m i c Flow Imaging Several important developments were necessary to arrive at the following proposed method for dynamic flow imaging. These include velocity phase mapping, the Oif-scan for inter-leaving images with dissimilar ghosting patterns and the above DIR technique.- In fact, the proposed dynamic flow imaging technique is a combination of the above techniques. The imaging sequence required combines the ideas of velocity field mapping with those of the OK-sc&n. First, two interleaved images are acquired using the OK-sc&n idea so that readout gradients, Gx and Gy are applied simultaneously, where Gy is applied with opposite polarity for each image as shown in Fig.4.2. Furthermore, the imaging sequences used are such that one image is velocity sensitized while the other image is velocity compensated. This is accomplished by applying appropriate gradients in the different directions. For simplicity, the velocity sensitization will be made to occur only along a single direction and hence the flow being imaged will be that along a single direction as well. The image processing involved for dynamic flow imaging combines the ideas of DIR as well as the assumption that the O if-scan can provide flow compensated and flow sensitized images, IComP{x,y) and Isens{x,y), with ghosting components, gComP(x,y) and gsens{x,y) re-105 spectively, that do not overlap at all. This is achieved by reversing the polarity of the Gy readout gradient in both Off-scan images so as to redirect the ghosting in different direc-tions. Once the two interleaved images are acquired, a subtraction: Isens(x,y) — Icomp{x,y) can be used to separate the time-averaged signal, Io(x,y) from gSens(x,y) a n d gComp{x,y), resulting in a combined ghost mask, gcs{x,y). This ghost mask will contain the required information from the ghosting in the compensated and sensitized images simultaneously but in different locations as a result of the Off-scan implementation. Therefore, the infor-mation contained in either gsens(x, y) or gCOmp(x, y) c a n be accessed by going to the correct location in gcs(x,y)- This would not be possible if the interleaved images have ghosts that are completely superimposed as those that were originally used for the 2pGEM method. For such images, the ghosting components that are overlaid on stationary signal are not as easily accessible. After the combined ghost mask is obtained, Eq.4.4 is applied to each pixel within a region-of-interest corresponding to the source of the flow and which will be referred to as ROI(Io). For now, this region can be approximated by a circle (since this is a good approximation for the cross-section of blood vessels) which is manually specified in terms of the position of the centre and its radius. However, any better ROI detection method could easily be implemented in the future. Each of the ghost masks yields a dynamic image at every pixel within ROI(Io), indicated by (x0,y0). The resulting dynamic images will be velocity sensitized, psens(xo,yo',t), a n d velocity compensated, pCOmp{xo, Vo] t), given by: PSenS{xo, yo] t) = ISens{xo, Vo) + 9cs(xo + nAxs, yQ + nAys)exp(inuj0t) n PcomP{xo, yo]t) = . IComp{xo, Vo) ~ 9cs{xo + nAxc, y0 + nAyc)exp(inu}0t} (4.6) n where n is the ghost order (from — nmax to +nmax) and U)Q can be calculated by Eq.4.5. As described for DIR, P 0 is the ghost spacing, given in pixels, along the phase encode direction. For this case of Off-scan images, the phase encode direction and the ghosting direction do not coincide but the P0 required to calculate UIQ still depends only on the spacing along the phase encode direction. This is because the interview motion, giving rise to the ghosting, modulates the signal between readout trajectories in &-space. Since these are related through a shift along the phase encode direction only, this is the direction of modulation and hence the shift in this direction reflects the periodicity of the motion. It is important to note that the time-averaged signal, fo, is identical for interleaved velocity sensitized and velocity compensated images only in the regions of stationary signal. This was discussed in section 3.4.4 where 2pGEM caused an unwanted signal cancellation of IQ in regions of the source of flow since f0 was not identical in the two interleaved images, in those regions. Therefore, since f0 in the ROI(IQ) is different for each image, when using 106 Eq.4.4, the appropriate value should be used for the reconstruction of each individual movie as indicated in Eq.4.6. Also, the incremental displacement step allowing to for access to each of the ghost orders in each interleaved image, appears along x and y directions rather than simply along the y direction as was the case in Eq.4.4. This is due to the obliquity of the ghosting distributions. More specifically, the displacement steps follow the oblique ghosting directions given by the angles as and ac for flow sensitized and compensated images respectively. These angles are defined with respect to the +a>axis as shown in Fig.3.15 and determined by the Oif-scan readout gradients. The displacement steps for the flow compensated and flow sensitized ghosting can thus be calculated in terms of the ghost spacing along the oblique direction, P0u, as: A x c A y c and A x s where P0 is the ghost spacing along the phase encode direction (assumed to be y here) and where the x displacement steps are related by: Axc = — Axs due to the relation between ac and a s . The two movies given by Eq.4.6 are complex, with real and imaginary components calculated separately. The first two terms of Eq.4.3 can be used to describe the contribution of the motion to the phase of the movies at a given time t. The phase of each of the movies as a function of time, 4>Sens(xo-> yo] t) and <t>Comp(xo, yo] t), can therefore be written: RTE RTE &sens(xo, yo] t) = <f>o(x0, y0] t) + jr0{x0, y0) / G s e n s ( t ' ) d t ' + j v 0 ( x 0 , y0) I t ' G s e n s ( t ' ) d t ' Jo Jo TE RTE G c o m p ( t ' ) d t ' + j v 0 ( x 0 , y 0 ) / t'Gcomp(t')dt(4.9) Jo where (f>o(xo, yo] t) represents the phase of the signal at time, t, induced by other factors than the motion and independent of the gradient waveforms. This term will be identical in both movies. The only difference between the phases in the two movies comes about because of the different gradient waveforms represented by the terms: G s e n s ( t ) and G c o m p ( t ) . If bipolar gradients are used to produce the two images as shown in Fig.4.2, the second term in Eq.4.9 will vanish for both sensitized and compensated movies but the third term will only vanish for the compensated movie, leaving the velocity dependence in the phase of the sensitized movie. Taking the phase difference, A(f)(x0,yo]t) = 4>senS(xo,yo]t) - (/>comP(xo,yo]t), will cause cancellation of the phase term, 4>(xo,yo]t), leaving only the velocity term: RTE A<P(x0,yo]t) = vo{x0,y0) / t ' G s e n s ( t ' ) d t ' . (4.10) Pobi • cos(ac) Pobi • sin(ac) = PQ P0bi • cos(as) Pobi • sin(as) Po-(4.7) (4.8) ycomP (x0, y0] t) = (po(xo, yo] t) + 7^0(^0, yo) Jo 107 A(f)(xo, yo', t) is thus proportional to the velocity of the spins at position (x0, y0) as a function of time and the sensitivity of the phase is controlled by the integral term which can be varied by varying Gsens(t). Since this flow sensitizing pulse will in general be given by a bipolar gradient waveform, the amplitude and duration of the bipolar lobes can be used to control the velocity sensitivity of a particular sequence. Since the flow is defined as being the integral of this localized velocity: v(xo, Vo) over the region ROI(I0), a summation of A(f)(x0, y0; t) over ROI(Io), should yield a dynamic flow profile, flow(t). For the purposes of this preliminary study, the relative shape of this profile is all that will be considered. 4.3.2 Determination of the Ghost Spacing As indicated in Eq.4.6, the fundamental frequency of the motion, o;0, must be known in order to obtain the required movies for dynamic flow imaging. Assuming this is not known a priori (as is usually the case), Eq.4.5 indicates that it can be determined from the spacing between ghosts along the phase encode direction, given by P0, and the total scan time, T 0, which is known for any given imaging sequence. Therefore, the determination of Po is an important step for dynamic flow imaging. Furthermore, Eq.4.6 also requires that the displacement steps, given by Eqs.4.7 & 4.8, be determined and these in turn depend on P 0 M , the oblique ghost spacing. Since PQu and P 0 are related through the known ghosting angles, ac or as, as given in Eqs.4.7 & 4.8, the determination of PQu is sufficient to determine the displacement steps and the fundamental frequency, given either ghost angle. This section describes an automated method that can be used to extract the ghost spacing, P0u from the images. In particular, it will be described for images with ghosting along an oblique direction but it can be applied, equally well, to the case of ghosting along the x- or y-axis. The only required input is the direction of the ghosting that can be entered as an angle, a, defined with respect to the -f-x-axis as in Fig.3.15. This angle is known for a given imaging sequence since it depends on the relation between the readout gradients: Gx and Gy which determine the fc-space trajectories during any readout period (see Chapter 1) and hence the corresponding ghosting direction (perpendicular to these trajectories). The profile, of the autocorrelation of a ghost mask image, along the direction of the ghosting and passing through the centre of the image, provides a very useful function from which the ghost spacing can usually be extracted quite easily. Since the autocorrelation of an image is the result of a convolution of the image with itself, any regular repetition of structures in an image may yield regularly spaced, off-centered peaks along a profile in the direction of the repetition. The combined ghost mask described above, gcs, will be made up of ghosting along two different directions. Its autocorrelation will therefore provide a 108 useful profile function along either one of these known directions. The off-centered peaks of either profile can be detected using a derivative operator and their spacing will correspond to the sought ghost spacing. This is the basis for the method of ghost spacing determination proposed here and given by the following steps. First, the combined ghost mask image, gcs, is created by a complex subtraction of the two interleaved images: ICOmp and Isens- The autocorrelation image, gcs*9cs, is then obtained by using the fact that a convolution (*) in image space corresponds to a multiplication in fc-space (Fourier transform space): gcs * gcs = FT-\FT{gcs) • FT{gcs)) (4.11) where FT and FT~l represent Fourier transform and inverse Fourier transform operators respectively. Then, one of the two ghosting directions is chosen: acornp or asens. For sim-plicity, the chosen direction will now be referred to simply as a. The central profile of the magnitude of gcs * gcs is taken along the chosen ghosting direction, a. The result is an oscillating function, p(r), with a relatively large central peak (at r = 0 which corresponds to the magnitude of gcs * gcs at the pixel given by (N/2, N/2), for an N x N image matrix). Every step in r corresponds to a step of r • cos(a) and r • sin(a) in the x and y direc-tions of the image respectively. If the ghosts are distinctly distributed, as is sometimes the case for flow-related ghosted images (as opposed to respiratory ghosted images), secondary peaks will be apparent in the associated p(r). Examples of a flow-related ghosted image, its autocorrelation and the corresponding profile are shown in Fig.4.3. A derivative operator, applied on the profile, p(r), can be used to produce another curve: pr(r) = dp(r)/dr. This curve, pr(r), is examined from the center outward in both directions. The first time the derivative vanishes, on either side of the central and overall maximum, corresponds to the first minimum. The second time the derivative vanishes, in a given direction away from the center, corresponds to the first secondary maximum in that direction. More generally, the 2nth time dp(r)/dr — 0 will correspond to the nth secondary maximum in that direction, where n is a positive integer. The ghost spacing will then be given by the spacing between two maxima, given by the difference between the corresponding r values: r; and ri+\ for the ith and (i + l)th maxima. There are several positions along pr (r) which can be used to determine the ghost spacing so an average of the several results is possible in order to improve the estimate. This ghost spacing determination technique was performed on the image shown in Fig.4.3 yielding P0u = 20.2 which correlates well with what is observed in the images. This value can be further checked since it should predict the period of the motion, T, which, for this experiment, was known since it was determined by the pump: T = 840 ms. P0bi was hence 109 Figure 4.3: Figure showing a magnitude ghosted image and the associated autocorrelation function with associated profiles, (a) shows the magnitude of one of the images, resulting from the OK-scm. (b) shows the magnitude of the ghost mask, gcs, resulting from the complex difference between the interleaved, ghosted, flow sensitized and compensated images. For these images, phase encoding is along the x direction and the frequency encoding is along the y direction hence the gradients are related by: Gy = 2GX, giving a = 26.6°. (c) is the autocorrelation of the image shown in (b). (d) is a plot of p(r), the central profile of image (c) along a = 26.6°. The actual profile has been smoothed a few times using a median filter in order to eliminate noise-induced spikes, (e) is pr(r), the derivative of the profile p(r) shown in (d). 110 used to determine T in the following way. For this case, the phase encode was along the x direction and the frequency encode was along the y direction so the gradients were related by: Gy = 2GX. This gives a = 26.6° for the image shown in Fig.4.3 and the cosine operation is used to relate P 0 to Pobi giving: P 0 = 20.2 • cos(26.6°) = 18.06. Eq.4.5 and the fact that LO0 = 2ir/T gives: T = T o / P 0 = (15.36/18.06) s = 850 ms. The ghost spacing determination therefore yielded the expected result within expected experimental error (< 2%). This technique has also been successfully applied on several other ghosted images. How-ever, in practice, it was found that applying a median filter in order to smooth the p(r) profiles was necessary in order to eliminate small noise-induced local extrema. 4.3.3 Experimental Set-up The goal was to qualitatively compare the resulting dynamic flow profiles obtained using the proposed fast method to a well-established clinical standard. In this case, VENE on. a 1.5 T scanner (EDGE, Picker International, Cleveland, OH, USA) was used, whereby prospective cardiac gating or triggering is used to collect interleaved flow sensitized and flow compensated images at several time frames within a given cycle. The resulting time series of flow sensitized and compensated images were then processed using a calibrated, semi-automated software package (AFLOW) in order to extract the wanted dynamic flow profile, defined as the changes in the integrated velocity, over a specified region-of-interest (ROI), as a function of time. Data was collected for a controlled and simplified set-up as well as for the realistic case of a human volunteer. The controlled set-up involved a programmable pump (UHDC flow system, Quest Imaging, London, Ontario, Canada) which was used to drive water through a hose in a determined, periodic way. The periodicity is quite accurate and an electric signal is generated by the pump, mimicking the electrical signal generated by the heart (ECG) to allow for the implementation of cardiac gated imaging. The flow pattern created inside the hose is initially unknown although the input at the pump is controlled. This is why a standard method of obtaining the actual dynamic flow profile should be used simultaneously for comparison with the result obtained using the fast method proposed. The pulse sequence used to collect the O if-scan data was created by modifying already existing interleaved flow compensated and flow sensitized gradient echo sequences. The flow sensitization of this sequence occurs along the z direction. The main modification involved the addition of a readout gradient, along the phase encode direction, of opposite polarities in both interleaved images. Unfortunately, this was accomplished on the chosen template sequence at the expense of flow compensation in the x and y directions. It is hoped that 111 at least for the controlled pump experiments, the flow is primarily in the z direction due to proper placement of the hose and hence this will not be a major source of contamination in the images. The resulting data was decoded using in-house software developed for restoring OK-scan images. However, this study is not dependent on this particular implementation. Similar experiments could be performed using any other available sequences resulting in dissimilar ghost patterns for the interleaved, flow compensated and sensitized images. Also, any other standard method of obtaining dynamic flow profiles could be used as a basis for comparison. The goal of this study, in its preliminary stage, is qualitative: the absolute value of the resulting dynamic flow profiles is not important. The attempt is to reproduce the relative shape of the profile obtained by a clinical standard, with the proposed fast OK-scan dynamic flow imaging. A quantitative comparison of the profile magnitudes and mean flow in time, would be the next step of the study. 4.4 Practical Considerations The complex form for flow sensitized and flow compensated movies can be created using the ideas of DIR as described by Eq.4.6. However, in order to progress from DIR to dynamic flow imaging two important steps must be taken. First, phase information must be extracted from the movies. In particular, it is the phase difference between the two flow movies that is of importance. Also, in order to obtain a dynamic flow profile, an integration of the flow throughout the ROI(I0) is necessary. Neither of these steps are trivial since the operations involve phase which behaves very differently from magnitude [61]. The most important distinction between the behaviour of phase and that of magnitude is due to the fact that the phase is only unambiguously defined within a closed interval of [—7r,7r) as discussed in Chapter 2. Problems arising due to this ambiguity in phase are well appreciated in the work presented in Chapter 2 and are of relevance again here. In practice, the phase difference between two complex numbers, c\ — Mie101 and c2 = M2e* ('2, can be calculated directly, without first calculating the phase for each individually. This is accomplished by multiplying c\ by the complex conjugate of c2, written c2: c, • c* = M1M2ei^B (4.12) since the phase difference, A9 = 8X — #2, is given by the phase of the above expression. Given any two complex numbers in terms of their real and imaginary parts: cx = R\ + iJx and c2 = i? 2 + the above multiplication can be performed by: d • c* = [(i?XJR2) + [Jrh)} + iKJM - (RM] (4.13) 112 where the real and imaginary components of the result are given between square brackets. From this, the phase of C\ -c2, equal to A9, can be calculated using an arctangent operation: Unfortunately, if the phase difference is determined this way, an implicit assumption is being made: the phase difference is assumed to be contained within [—7r,7r). This results since there are always two possible angles that describe the position of one vector relative to another and the arctangent operation results in the principal phase value (see Chapter 2) which is the smallest of the two angles. This is further complicated by the fact that any multiple of 2-7T can always be added or subtracted depending on the number of wrap arounds. There are therefore at least two different ways of determining the phase difference be-tween two complex numbers. Either indirectly, by first calculating the principal phase value (contained within ••[—7r, 7r)) for each individually and then subtracting the two results or directly, by using the complex conjugate multiplication described above. The difference between these two methods is in the allowable range of phase difference. Due to the control of velocity sensitization given by the strength and duration of the bipolar flow sensitizing pulse as given in Eq.4.10, it is reasonable to assume that no wrap arounds occur between the phase of the flow sensitized and compensated movies. This argument could be extended to assume that the phase difference due to the constant velocity of the object being imaged, should be contained within the interval [—7r,7r). It is therefore reasonable to use the complex conjugate multiplication described above to extract A8 at each pixel. This is the implementation used in the following. The desired dynamic flow profile is then calculated by simply integrating these phase differences over the ROI(I0). The process of phase integration merits further consideration. Phase is very sensitive to the effects of noise. Since noise is, in general, of smaller magnitude than relevant signal, it will not have such an effect on the magnitude of the resulting signal. However, low signal pixels associated with noise are more likely to have any phase associated with them and this can greatly alter the result of the integration. This is another important consideration when choosing the exact implementation for the proposed dynamic flow imaging since the result may be more or less sensitive to noise depending on the order of the steps taken. For example, if the real and imaginary contributions are first integrated and the resulting phase is then taken, the result would not be the same. This order of operations would present the same benefits of averaging noise that occurs in magnitude calculations since relative attenuation of noise would result due to the relative magnitudes of noise and relevant signal. However, this order of obtaining an average phase is not equivalent to first taking the phase (4.14) 113 at each pixel and then integrating the phase over all the pixels in ROI(Io)- An appreciation of the sensitivity of phase to such operations and their order has been addressed by other researchers as well [56, 61]. Ultimately, since it is an integration of the velocity at each pixel that is wanted, it seems more reasonable to integrate the individual phases since the flow at each pixel is given by the phase at each pixel. This is the rationale that was used and therefore that is the order of steps used for phase integration that was implemented in the following. However, it is important to note that the result will be very sensitive to noise and therefore the ROI(Io) determination is important. For now, this is approximated as a circle of prescribed radius and location and this may be the source of some error. In conclusion, the effects of noise must be further explored. Another consideration of the implementation is that the method relies on the fact that ghosting in Icomv and Isens are completely not overlapping with each other such that all the necessary information, from gcomp and g3ens, is readily accessible in gcs. However, this will depend on the number of ghosted regions, their relative positions and the exact ghosting direction. Fig.4.4 demonstrates two cases where the ghosts of the interleaved images overlap significantly, possibly affecting the implementation of the proposed method. Luckily, the set of interleaved images collected in this study were obtained twice, with swapped phase and readout directions. This resulted, in general, in at least one set of images with nicely distributed ghosts for each experimental set-up. The last consideration that will be addressed is the problem of defining absolute time. The proposed method uses an analytical expression, given by Eq.4.6, to incorporate the time-dependence. This expression incorporates information taken directly from the image to describe the time-averaged, DC, component of the modulation as well as the harmonics. The production of these harmonics (ghosts) and the DC component (I0) can be explained by looking at fc-space. Therefore, the time origin t = 0, generated by Eq.4.6 is defined as the instance when the central A;-space data is acquired. Without the use of the ECG signal, it is not clear as to which instant this corresponds to in the cardiac cycle. In contrast, this absolute time definition is not ambiguous for the movie created by several snap shots in time, as is done for cardiac gated images. This is because in these cases, time is defined in terms of a specific trigger from the ECG signal. It is therefore important to note that although a similar shape to the resulting dynamic profile is hoped to be achieved by the proposed technique and cardiac gating, a shift along the time axis may be necessary in order to observe a match. This should be allowed due to this ambiguity in the definition of absolute time. There are obvious considerations that must be further investigated in order to best 114 Figure 4.4: Figure showing two examples of a set OK-scan interleaved flow sensitized and flow compen-sated images, with ghosts that interfere. One set-up consisted on two stationary water bottles and a looped hose with water being pumped at a controlled rate, (a) and (b) show the magnitude interleaved images and (c) shows the magnitude of their complex difference. The other set of images, shown in (d) and (e) correspond to the cross-section of the abdomen of a human volunteer. Two major vessels can be seen in close proximity and their respective ghosts are seen to interfere with one another. 115 extract the dynamic flow profile. The purpose of this section was to describe and explain the rationale for the choice of implementation used to produce the preliminary results. In doing so, the complexity of the apparently simple idea has hopefully become apparent. The source of the complexity lies in the fact that the velocity is encoded in the phase of the signal and phase is a very ambiguous and sensitive measure. 4.5 Preliminary Results Due to all the considerations mentioned in section 4.4, it must be emphasized that the results presented in this section are only preliminary. The exact implementation of the basic idea for dynamic flow imaging has many unexplored options such as varying the many parameters involved in the pulse sequence to best optimize the distribution of the ghosts and varying the processing that averages the time-varying velocity within a given ROI. This means that many things have not been optimized for the. following results and yet the dynamic flow profiles resulting from some of the processed images are quite promising. These results are therefore significant since they can be used to justify further research in this area. The most promising results obtained so far are for the controlled experimental set-up with pump and hose. The volunteer images resulted in too much interference between the ghosts as shown in Fig.4.4. Two cases will be shown as representative of the good results observed. Both cases are based on the same experimental set-up consisting of a hose with pumped water and some stationary bottled water. In both cases, the hose is looped around so that the cross-sectional image contains two views of the hose, one with water flowing in the positive z direction and the other with it flowing in the opposite direction. This means that two ROI(IQ)S are present in the images and each one has its associated ghosts. While this did make it more of a challenge to obtain the two interleaved OK-sc&n images with ghosts that do not interfere significantly, it also served as a very useful check. Although the main goal is to qualitatively reproduce the profile obtained using VENE and AFLOW, by setting up the hose in this looped form, another important check was made available. This check compares the dynamic flow profiles observed in each of the two cross-sections of the hose. Since the flow profiles can be expected to be very similar in shape, only reversed in direction, a comparison of the resulting dynamic flow profiles should reflect this. Fig.4.5 shows the dynamic flow profiles obtained by performing the proposed processing on the shown set of OK-sc&n images. Dynamic flow profiles given by processing of the standard VENE sequence images were also obtained and are shown for comparison. The profiles behave qualitatively as expected: the profiles from the two tube cross-sections are 116 similarly shaped but reversed and they are in qualitative agreement with the shapes obtained by the standard VENE sequence. In Fig.4.6, the proposed processing was also shown to give the expected, reversed profiles, for the two ROI(I0)s even though there was some ghost interference for orders greater than three as demonstrated in Fig.4.4. Unfortunately, at the time of these experiments, a full cycle was not simultaneously acquired using the standard VENE sequence. However, since the pump and hose set-up was identical to the case shown in Fig.4.5, similar flow profiles can be expected so the profiles shown in (d) and (e) of that figure can be used as a basis for comparison with those shown in Fig.4.6 (c). These are, in fact, qualitatively similar. 4.6 Discussion Sz Conclusion Although the dynamic flow imaging described above shows some promising, preliminary results, it relies on the Off-scan type of data acquisition (producing dissimilar ghosting patterns) which is not readily available. Another possibility exists, for a simpler imple-mentation of the data acquisition. This method would not rely on the Off-scan, but its robustness is more uncertain since it relies more on processing and in particular, on the 2pGEM method of ghost suppression. Although this idea has not been fully explored, it will be discussed here. The idea is to acquire the two interleaved images, f c o m p and I sens, in a standard way, with a readout gradient exclusively along the x direction. In such a case, the ghosts will not be obliquely directed, but rather along the y direction. If there is sufficient difference between the ghosts of the two images due to the flow compensation and sensitization, then the conditional 2pGEM method can be applied to extract f0. Otherwise a Ak can be used as before. Removing the fo component from the original images should result in the ghost masks: gCOmV and gsens respectively. These ghost masks will then provide the necessary information to create flow compensated and sensitized movies once the ghost spacing is known. Either one of these ghost masks can be used to determine the ghost spacing as was described above, given the appropriate ghosting direction. Therefore, such a method would not require the Off-scan type of implementation and the same dynamic flow information could be obtained. However, this method has some potential problems. First of all, it relies on the assump-tion that the result of the 2pGEM method exactly represents fo. This is used to remove the f0 component from the originally interleaved, ghosted images in order to extract the ghost masks. Since the 2pGEM method actually produces a result that best minimizes the gradient energy criterion, it is not an absolute extraction of the f0. In fact, as has been 117 Figure 4.5: Figure comparing dynamic flow profiles obtained by the proposed technique and using a standard clinical method. The experimental set-up consists of a stationary water bottle and a looped hose with water being pumped at a controlled rate. The magnitudes of the interleaved, OK-scan flow sensitized and flow compensated images are shown in (a) and (b). The two ROI(Io)s, corresponding to the two hose cross-sections, can be distinguished since they appear brighter than their ghosts. One hose cross-section contains water flowing in the +z direction while the other contains water flowing in the — z direction. The profiles resulting from the proposed method, for both ROI(Io)s, are shown in (c) while those resulting from the standard V E N E sequence are shown in (d) and (e). (f) shows a match between the two different results obtained for the same ROI: the corresponding profile from (c) has been time shifted, scaled and overlaid on the profile in (d). 118 Figure 4.6: Figure comparing dynamic flow profiles obtained by the proposed technique in two cross-sections of the same hose. The experimental set-up consists of two stationary water bottles and a looped hose with water being pumped at a controlled rate. The magnitudes of the interleaved, Off-scan flow sensitized and flow compensated images are shown in (a) and (b). The two ROIs, corresponding to the two hose cross-sections, can be distinguished since they appear brighter than their ghosts. One hose cross-section contains water flowing in the +z direction while the other contains water flowing in the —z direction. The profiles resulting from the proposed method, for both ROIs, are shown in (c). discussed in the previous chapter, the 2pGEM method provides a balance between noise amplification and ghost suppression. This means that in some regions, specially those with low signal, the result is more dominated by noise which drives the result closer to a simple average of the two images, rather than I0. Since the higher order harmonics of motion tend to result in ghosts of smaller signal intensity, the result of the 2pGEM in those regions is not expected to represent I0. The ghost masks produced by removal of the result of the 2pGEM will thus not generally be a good representation of the higher order harmonics, and if used to obtain the dynamic flow, inaccurate results would be expected. The ideas discussed in this section need to be further explored. In particular, the prac-tical considerations associated with the implementation of the technique need to be studied in order to obtain a more optimized implementation. The hope is, nonetheless, that the studies presented in this chapter demonstrate the feasibility of applying the DIR to velocity phase encoded images and that the promising preliminary results can be used to motivate further development of the technique. Once this technique is better established, the 2pGEM method of ghost suppression may prove very useful in extending its applicability. 119 Chapter 5 Conclusion This thesis consisted of three distinct projects sharing the common focus of MRI artifacts. In particular, the phase wrap artifact and the motional ghosting artifact, which both greatly limit the diagnostic potential of MRI, were studied. The goal of the work was to better understand and characterise these artifacts in order to devise methods of eliminating or reducing their effects on the images. The knowledge that resulted from their study, as well as the tools developed for dealing with them, allowed for the realisation of a fast dynamic flow imaging technique where the phase information contained within the ghosting artifact is extracted. The first project led to an improved understanding of phase maps by a better character-isation of the pole structures. Although the aim here was to exploit this knowledge in order to successfully phase unwrap, such an understanding of phase maps could be extended to many other applications. Since poles are always present in the noise (in fact, there are one third of the number of noise pixels) their localization could be used to distinguish between noise and signal in other applications such as tissue segmentation or statistical estimates of signal-to-noise ratios. Furthermore, a quality map was developed which makes use of the connectivity of the dipole structures in phase maps. This map is quite unique in that it highlights the relative phase gradients contained in phase maps only in areas that are affected by poles. This is believed to be useful information about a phase map since the poles are the origin of the phase unwrap problem and possibly other phase-related problems. Being able to obtain a quality map focusing on those regions proved very useful as a guide for phase unwrapping and it is likely that it could be useful in other applications. The method of phase unwrapping that resulted from this study was shown to be quite successful in unwrapping phase maps for which conventional methods would have failed. It was also shown to be quite robust to the effects of noise. The main limitation of the method was due to its two-dimensional nature and an extension to three dimensions was proposed. The second project led to several improvements of an already existing ghost suppression 120 technique which relies on the input of two ghosted images produced from a very specific data acquisition scheme. It then combines the two input images with parameters selected according to the criterion of minimization of the gradient energy of the resulting image. The development of an analytical form for gradient energy minimization (GEM) was a great stepping stone in this research. Not only did it result in a far better implementation of an already existing technique that became referred to here as lpGEM, but the simplicity, accuracy and efficiency of the minimization enabled the development of a more general form of ghost suppression. This method, called 2pGEM, lessened the restrictions of the two input images, allowing for a much broader application. The 2pGEM method, in turn, motivated the development of a new data acquisition scheme, called the OK-scan, where the ghosting patterns of two interleaved images are produced along different, oblique directions. The new data acquisition scheme allowed for the realisation of a third project, fast dynamic flow imaging, which is discussed below. Furthermore, the extended use of the 2pGEM method, which relies on the gradient energy minimization as a criterion to balance ghost suppression and noise amplification, served to support the use of the gradient energy as an indicative measure of image quality. Familiarity with the gradient energy of an image (or region) also resulted in two specific applications: general lpGEM for magnitude-only images and the conditional 2pGEM which depends on gradient energy relations. These will hopefully motivate further exploitation of this measure. Another product of the 2pGEM method was the ghost suppression method using one and a half or fewer excitations. This method was theoretically presented using simulations to demonstrate its ability. Although some limitations must still be overcome for it to become more practical, it may lead the way to other exciting advances in MRI by reducing the amount of necessary data collection. This hope is based on the following. The 2pGEM method of ghost suppression has been applied exclusively to motional ghost-ing artifacts in this thesis. Although this artifact is the focus of much of this thesis, another ghosting artifact was uncovered in the process: the undersampling ghost. In section 3.6, the undersampling ghost of the time-averaged part of the images, I0, is treated as quite distinct from the motional ghosting artifact. It is eliminated by simple addition of two complementary images, such as one resulting from sampling only the even lines of A;-space and the other from sampling only the odd lines. However, the 2pGEM method should be able to reduce this type of ghosting as effectively as the motional ghosting since the only restrictions it imposes on the images is that there be some significant difference (in phase and/or magnitude) between the ghosts. In fact, there is an inconsistency in section 3.6, between the treatment of the undersampling ghost of I0 and the undersampling ghosts of the motional ghosts. The latter appear in the intermediate images as input for the 2pGEM 121 method and are hence treated like the original motional ghosts. This emphasizes the fact that the 2pGEM method could be used to suppress undersampling ghosts. Similarly, this work could help benefit echo-planar images (EPI) contaminated by the N/2 ghost [54]. The reason this was not implemented nor presented in section 3.6, is due to the much discussed window size problem. The characteristic size of I0 is usually very different from the motional ghosts and, hence, the best possible results from the 2 p G E M method would require different implementations with varying window sizes and possibly shapes. Further-more, the ghost overlap problem is much more of an issue for undersampling ghosts of I0 than for motional ghosts due also to their larger coverage of the image space. Both these issues restricted the use of the 2pGEM method to the motional ghosts and their under-sampling counterparts. However, if the ghost overlap problem is overcome, suppression of undersampling ghosts by the 2pGEM method could become more practical. This would have the potential of greatly increasing the efficiency of data collection in MRI. As a result, the described method of motional ghost suppression could be further extended to less than a single excitation. The last project presented in this thesis, called fast dynamic flow imaging, involved the use of the motional ghosting artifact in velocity phase-encoded MR images, to reproduce the time-evolution of the integrated velocity of spins in a given ROI of an image. This project was developed as a way of combining several already existing methods such as dynamic image reconstruction (DIR) and velocity phase-encoded MRI. Its motivation was a direct result of the first two studies which focussed on phase maps as well as the motional ghosting artifact. A fast dynamic flow imaging method therefore appeared as a candidate, that could benefit immensely from the newly developed phase unwrapping method as well as the understanding of the information contained within the motional ghosting artifact. However, it was not until the realisation of the new data acquisition scheme, the OK-sc&n, that the dynamic flow imaging project could be undertaken. This is because the OK-scan provided a means of acquiring two interleaved images with dissimilar ghosting patterns and hence the potential of obtaining a useful ghost mask by simple subtraction. Although the preliminary results are promising, there remain many unanswered questions which will hopefully help guide future research. Finally, each of the techniques presented in this thesis could benefit from the input of more information. The phase unwrapping technique is limited in two-dimensions as discussed and could greatly benefit from incorporation of another dimension of information. The limitations of the 2pGEM method of ghost suppression have to do with ghost overlap and window size and shape selection. These limitations are local in that only certain regions of the images are affected by them. A third acquisition could help resolve this by providing 122 complementary information which could be used to help correct the affected regions. A similar argument could be used for the fast dynamic flow imaging technique which is also limited by the ghost overlap problem. A direction of future research in all these projects is thus to determine the best way to incorporate new information in order to overcome some of the limitations of the techniques in their present form. Ultimately, the new phase unwrapping method could be incorporated in the processing of the fast dynamic flow imaging. Also, another potentially better implementation of the fast dynamic flow imaging, which relies on the 2pGEM method of ghost suppression and is described in section 4.6, could be used. Incorporation of the new phase unwrapping technique and the 2pGEM method in the implementation of the fast dynamic flow imaging would tie together all three projects undertaken in this thesis. The work of this thesis resulted in some useful novel techniques for processing MRI data with artifacts. This was achieved by making extensive use of the phase contained within the MR images. Such a thorough investigation has also led to a better appreciation of the complexity of dealing with phase information. It is believed that the importance of the information encoded in the phase of MR images, as well as its potential use, are not fully appreciated. Hopefully, this thesis has emphasized the incredible potential of MRI due to its unique ability of encoding important information in the phase of the signal. It is also hoped that the development of the ideas throughout the thesis will aid in future exploitation of the phase and artifacts of an MR image. 123 Glossary of Abbreviations lpGEM : one-parameter gradient energy minimization 2pGEM : two-parameter gradient energy minimization 2D : two-dimensional 3D : three-dimensional 2DFT : two-dimensional Fourier transform CT : computerized tomography DC : direct current DIR : dynamic image reconstruction Eg : gradient energy ECG : electro-cardiogram FID : free-induction decay signal FGV : field-of-view FT : Fourier transform FT'1 : inverse Fourier transform Gx : magnetic field gradient along the x direction Gy : magnetic field gradient along the y direction Gz : magnetic field gradient along the z direction GEM : gradient energy minimization h : time-averaged signal TJO • weighted sum ISMRM : International Society of Magnetic Resonance in Medicine NMR : nuclear magnetic resonance NSA : number of signal averages MR : magnetic resonance MRI : magnetic resonance imaging OK-scan : oblique /c-space scan PSF : point-spread function RF : radio-frequency pulse ROI : region-of interest SNR : signal-to-noise ratio TR : repetition time TE : echo time 124 Bibliography [1] P.C. Lauterbur. Image formation by induced local interactions: examples employing nuclear magnetic resonance. Nature, 242:190-191, 1973. [2] I.I. Rabi, J.R. Zacharias, S. Millman, and P. Kusch. A new method of measuring nuclear magnetic moment (letter). Physical Review, 53:318, 1938. [3] I.I. Rabi, S. Millman, P. Kusch, and J.R. Zacharias. The molecular beam resonance method for measuring nuclear magnetic moments. Physical Review, 55:526-535, 1939. [4] Q.-S. Xiang. Introduction to MRI. Class notes: UBC-PHYS404, 2000. [5] F. Bloch, W.W. Hansen, and M. Packard. Nuclear induction. Physical Review, 69:127, 1946. [6] F. Bloch. Nuclear induction. Physical Review, 70:460-474, 1946. [7] E.M. Purcell, H.C. Torrey, and R.V. Pound. Resonance absorption by nuclear magnetic moments in a solid. Physicsl Review, 69:37-38, 1946. [8] M.L. Wood. Fourier imaging. In Magnetic Resonance Imaging, volume 1, chapter 2. Mosby Year Book, St.Louis, MO, 2nd edition, 1992. [9] E.L. Hahn. Spin echoes. Physical Review, 80:580-594, 1950. [10] R.M. Henkelman and M.J. Bronskill. Artifacts in magnetic resonance imaging. Reviews of Magnetic Resonance in Medicine, 2(1):1—126, 1987. [11] D.C. Ghiglia and M.D. Pritt. Two-Dimensional Phase Unwrapping. Wiley, New York, 1998. [12] S. Chavez and Q.-S. Xiang. 2D phase unwrapping based on dipole connections. In Pro-ceedings of the International Society for Magnetic Resonance in Medicine, 5th Scientific Meeting and Exhibition, page 2053. ISMRM, April 1997. 125 [13] S. Chavez, Q.-S. Xiang, and L. An. Understanding phase maps in MRI: A new cut-line phase unwrapping method. IEEE: Transactions on Medical Imaging, 2000. re-submitted in final form, August 2000. [14] S. Chavez and Q.-S. Xiang. Efficient ghost suppression by analytical gradient energy minimization. In Proceedings of the International Society for Magnetic Resonance in Medicine, 6th Scientific Meeting and Exhibition, page 132. ISMRM, April 1998. [15] S. Chavez and Q.-S. Xiang. Improved ghost suppression by two-parameter gradient en-ergy minimization. In Proceedings of the International Society for Magnetic Resonance in Medicine, 7th Scientific Meeting and Exhibition, page 1999. ISMRM, May 1999. [16] S. Chavez and Q.-S. Xiang. A quantitative evaluation of ghost suppression algorithms. In Proceedings of the International Society for Magnetic Resonance in Medicine, 8th Scientific Meeting and Exhibition, page 1691. ISMRM, April 2000. [17] S. Chavez and Q.-S. Xiang. Fast ghost suppression using 1.5 or fewer excitations. In Proceedings of the International Society for Magnetic Resonance in Medicine, 7th Scientific Meeting and Exhibition, page 239. ISMRM, May 1999. [18] LR. Young and G.M. Bydder. Phase imaging. In Magnetic Resonance Imaging, vol-ume 1, chapter 9. Mosby Year Book, St.Louis, MO, 2nd edition, 1992. [19] R. Kramer and O. Loffeld. A novel procedure for cutline detection. International Journal of Electronic Communications, 50(2):112-116, March 1996. [20] S. Nakadate and H. Saito. Fringe scanning speckle-pattern interferometry. Applied Optics, 24(14):2172-2180, 1985. [21] J.R. Buckland, J.M. Huntley, and S.R.E. Turner. Unwrapping noisy phase maps by use of a minimum-cost-matching algorithm. Applied Optics, 34(28):5100-5108, August 1995. [22] J.R. Reichenbach, R. Venkatesan, D.A. Yablonskiy, M.R. Thompson, S. Lai, and E.M. Haacke. Theory and application of static field inhomogeneity effects in gradient echo imaging. Journal of Magnetic Resonance Imaging, 7(2):266-279, March/April 1997. [23] T.R. Judge and P.J. Bryanston-Cross. A review of phase unwrapping techniques in fringe analysis. Optics and Lasers in Engineering, 21:199-239, 1994. [24] J. Strand and T. Taxt. Performance evaluation of two-dimensional phase unwrapping algorithms. Applied Optics, 38(20):4333-4344, July 1999. 126 [25] C. DeVeuster, P. Slangen, Y. Renotte, L. Berwart, and Y. Lion. Disk-growing algorithm for phase-map unwrapping: application to speckle interferograms. Applied Optics, 35(2):240-247, January 1996. [26] R.M. Goldstein, H.A. Zebker, and C.L. Werner. Satellite radar interferometry: two-dimensional phase unwrapping. Radio Science, 23(4)-.713-720, July/August 1988. [27] J.M. Huntley. Noise-immune phase unwrapping algorithm. Applied Optics- letters to the editor, 28(15), 1989. [28] D.J. Bone. Fourier fringe analysis: the two-dimensional phase unwrapping problem. Applied Optics, 30(25):3627-3632, September 1991. [29] D.P. Towers, T.R. Judge, and P.J. Bryanston-Cross. Automatic interferogram analysis techniques applied to quasi-heterodyne holography and ESPI. Optics and Lasers in Engineering, 14:239-281, 1991. [30] N.H. Ching, D. Rosenfeld, and M. Braun. Two-dimensional phase unwrapping using a minimum spanning tree algorithm. IEEE: Transactions on Image Processing, 1(3):355-365, July 1992. [31] M.Takeda and T. Abe. Phase unwrapping by a maximum cross-amplitude spanning tree algorithm: a comparative study. IEEE: Transactions on Image Processing, 1992. [32] Q.-S. Xiang. Temporal phase unwrapping for CINE velocity imaging. Journal of Magnetic Resonance Imaging, 5(5):529-534, September/October 1995. [33] L. An, Q.-S. Xiang, and S. Chavez. Fast minimum spanning tree method for phase unwrapping. IEEE: Transactions on Medical Imaging, 2000. [34] T.J. Flynn. Consistent 2-D phase unwrapping guided by a quality map. In Proceed-ings of the International Geoscience and Remote Sensing Symposium (IGARSS), pages 2057-2059. IEEE, New York 1996. [35] G. Scott and M.L.G. Joy. Guaranteed phase unwrapping without noise thresholding. In Proceedings of the Society for Magnetic Resonance in Medicine, 10th Annual Meeting and Exhibition, 10, page 747. SMRM, 1991. [36] Q.-S. Xiang, S. Chavez, K. Whittall, and G. Culham. ROI definition using score-guided-erosion. In Proceedings of the International Society for Magnetic Resonance in Medicine, 5th Scientific Meeting and Exhibition, page 2027. ISMRM, April 1997. 127 [37] G.L. Nayler, D.N. Firmin, and D.B. Longmore. Blood flow imaging by cine magnetic resonance. Journal of Computer Assisted Tomography, 10(5):715-722, 1986. [38] M.L. Wood and R.M. Henkelman. MR image artifacts from periodic motion. Medical Physics, 12(2): 143-151, 1985. [39] W.Ff. Perman, P.R. Moran, A. Moran, and M.A. Bernstein. Artifacts from pulsatile flow in mr imaging. Journal of Computer Assisted Tomography, 10(3):473-483, 1986. [40] E.M. Haacke and J.L. Patrick. Reducing motion artifacts in two-dimensional fourier transform imaging. Magnetic Resonance Imaging, 4:359-376, 1986. [41] M.L. Wood. Respiratory artifacts: Mechanism & control. In D.M. Grant and R.K. Harris, editors, Encyclopedia of Nuclear Magnetic Resonance, pages 4139-4146. John Wiley k Sons, 1996. [42] Q.-S. Xiang and R.M. Henkelman. Motion artifact reduction with three-point ghost phase cancellation. Journal of Magnetic Resonance Imaging, 1:633-642, 1991. [43] Q.-S. Xiang and R.M. Henkelman. Two-point interference method for suppression of ghost artifacts due to motion. Journal of Magnetic Resonance Imaging, 3:900-906, 1993. [44] Q.-S. Xiang and S. Chavez. Oblique k-space scan for motion artifact suppression. In Proceedings of the International Society for Magnetic Resonance in Medicine, 8th Scientific Meeting and Exhibition, page 398. ISMRM, April 2000. [45] Q.-S. Xiang and R.M. Henkelman. Weighted average and its application in ghost artifact reduction. 9th SMRI, page 222, 1991. [46] D.G. Kruger, G.S. Slavin, R. Muthupillai, R.C. Grimm, and S.J. Riederer. An orthogo-nal correlation algorithm for ghost reduction in MRI. Magnetic Resonance in Medicine, 38:678-686, 1997. [47] R.S. Hinks, J.A. Derbyshire, and P.L. LeRoux. Virtual frequency encoding: Achieving symmetry in K-space. In Proceedings of the International Society for Magnetic Res-onance in Medicine, 6th Scientific Meeting and Exhibition, page 570. ISMRM, May 1998. [48] Q.-S. Xiang. Reducing motion and off-resonance effects with half-cartesian k-space sam-pling. In Proceedings of the International Society for Magnetic Resonance in Medicine, 6th Scientific Meeting and Exhibition, page 243. ISMRM, May 1998. 128 [49] M.L. Wood. Ineffectiveness of averaging for reducing motion artifacts in half-fourier MR imaging. Journal of Magnetic Resonance Imaging, 1:593-600, 1991. [50] B.Madore and R.M. Henkelman. A new way of averaging with applications to mri. Medical Physics, 23(1):109-113, January 1996. [51] D. Atkinson, D.L.G. Hill, P.N.R. Stoyle, R E . Summers, and S.F. Keevil. Automatic correction of motion artifacts in magnetic resonance images using an entropy focus criterion. IEEE: Transactions on Medical Imaging, 16(6):903-910, 1997. ' [52] K.P. McGee, A. Manduca, J.P. Felmlee, S.J. Riederer, and R.L. Ehman. Image metric-based correction (autocorrection) of motion effects: Analysis of image metrics. Journal of Magnetic Resonance Imaging, 11:174-181, 2000. [53] P.M. Margosian, G. DeMeester, and Haiying Liu. Partial fourier acquisition in MRI. In D.M. Grant and R.K. Harris, editors, Encyclopedia of Nuclear Magnetic Resonance, pages 3462-3467. John Wiley & Sons, 1996. [54] M.H. Buonocore and L. Gao. Ghost artifact reduction for echo planar imaging using image phase correction. Magnetic Resonance in Medicine, 38:89-100, 1997. [55] W.G. Bradley. Flow phenomena. In Magnetic Resonance Imaging, volume 1, chapter 11. Mosby Year Book, St.Louis, MO, 2nd edition, 1992. [56] P.R. Moran. A flow velocity zeugmatographic interlace for NMR imaging in humans. Magnetic Resonance Imaging, 1:197-203, 1983. [57] D.J. Bryant, J.A. Payne, D.N. Firmin, and D.B. Longmore. Measurement of flow with NMR imaging using a gradient pulse and phase difference technique. Journal of Computer Assisted Tomography, 8(4):588-593, 1984. [58] P. van Dijk. Direct cardiac NMR imaging of heart wall and blood flow velocity. Journal of Computer Assisted Tomography, 8(3):429-436, 1984. [59] D.N. Firmin, G.L. Nayler, R.H. Klipstein, S.R. Underwood, R.S.O. Rees, and D.B. Longmore. In vivo validation of MR velocity imaging. Journal of Computer Assisted Tomography, 11 (5):751-756, 1987. [60] Q.-S. Xiang and R.M. Henkelman. Dynamic image reconstruction: MR movies from motion ghosts. Journal of Magnetic Resonance Imaging, 2:679-685, 1992. 129 [61] C A . Hamilton. Analysis and Improvements in Magnetic Resonance Phase Contrast Flow Imaging. PhD thesis, Department of Electrical and Computer Engineering: North Carolina State University, Raleigh, North Carolina, February 1992. 130 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085490/manifest

Comment

Related Items