Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Viscoelastic imaging methods using acoustic radiation force Frew, Samuel Mark 2010

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

Item Metadata


24-ubc_2011_spring_frew_samuel.pdf [ 5.71MB ]
JSON: 24-1.0071534.json
JSON-LD: 24-1.0071534-ld.json
RDF/XML (Pretty): 24-1.0071534-rdf.xml
RDF/JSON: 24-1.0071534-rdf.json
Turtle: 24-1.0071534-turtle.txt
N-Triples: 24-1.0071534-rdf-ntriples.txt
Original Record: 24-1.0071534-source.json
Full Text

Full Text

Viscoelastic ImagingMethods Using Acoustic Radiation Force by Samuel Mark Frew BE (Hons), University of Canterbury, 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2010 © Samuel Mark Frew 2010 Abstract Elastography is a method of imaging the viscoelastic properties of soft tissue and similar media. It can be used for identifying various forms of cancer and fibrosis in soft tissue. Elastography can thus aid in diagnosis of such diseases, and in treatment procedures, such a biopsy needle guidance. Acoustic radiation force (ARF) imaging is one form of elastography whereby the ARF of focused ultrasound displaces the medium at an internal, localised position. The same ultrasound system can be used to monitor the displacement, from which the viscoelastic properties can be recovered. This work presents a new ARF imaging method, called axial relaxation imag- ing, that uses transfer function techniques. The relaxation response at a point in the medium is monitored after a period of ARF application. By assuming the response corresponds to a negative step response of a linear system, the rela- tive force-displacement transfer function is computed. This is then related to a Voigt model over a range of frequencies to obtain a relative elastic parameter and a frequency-dependent viscous parameter governed by a power law. The method was applied to homogeneous phantoms made with different gelatine concentrations. Relative elastic parameters of 1, 2.3 and 4.4 and relative viscous parameters of 1, 1.8 and 3.0 were obtained for 2, 3 and 4 wt% gelatine concentrations, respectively, with consistent results between phantoms of the same type and agreement with values from other estimation techniques. The viscous power law frequency dependencies were governed by flow index values of –0.10, –0.14 and –0.18, respectively. The good separation between parameters in the results shows the method holds potential for application to tissue charac- terisation. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation and Overview . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Ultrasonic Imaging . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2 Fundamentals of ARF . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Review of ARF Imaging Techniques . . . . . . . . . . . . . . . . . . . 8 1.3.1 Static Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.2 Transient Methods . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.3 Periodic Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3.4 Perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 Ultrasound Platform and Probe . . . . . . . . . . . . . . . . . . . . . 15 2.2 Texo Sequence Generating Software . . . . . . . . . . . . . . . . . . . 17 2.2.1 Texo SDK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 iii Table of Contents 2.2.2 ARF Imaging Software . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.3 Sequence Parameters . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.4 Axial and Lateral Relaxation Imaging Sequences . . . . . . . 23 2.3 Data Processing Software . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3.1 Data Import . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3.2 Displacement Estimation . . . . . . . . . . . . . . . . . . . . . 27 2.3.3 Axial Relaxation Imaging: Transfer Function Estimation . . . 30 2.3.4 Axial Relaxation Imaging: Parameter Estimation . . . . . . . 32 2.3.5 Lateral Relaxation Imaging: Shear Wave Velocity Estimation 36 2.3.6 Lateral Relaxation Imaging: Parameter Estimation . . . . . . 38 2.3.7 Display Conversion . . . . . . . . . . . . . . . . . . . . . . . . 38 2.4 Transmission Medium . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.2 Phantom Construction . . . . . . . . . . . . . . . . . . . . . . 40 2.5 Development Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.5.1 Heating and ARF . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.5.2 Generating Pushing Pulses . . . . . . . . . . . . . . . . . . . . 44 3 Safety Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1.1 Heating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1.2 Cavitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.1.3 Acoustic Output Parameters . . . . . . . . . . . . . . . . . . . 48 3.2 Acoustic Output Experiment . . . . . . . . . . . . . . . . . . . . . . . 50 3.2.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3 Probe Heating and Protection . . . . . . . . . . . . . . . . . . . . . . 55 4 ARF Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 iv Table of Contents 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5 Homogeneous Phantom Experiments . . . . . . . . . . . . . . . . . . . . 68 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.3 Axial Relaxation Imaging Results . . . . . . . . . . . . . . . . . . . . . 69 5.3.1 Displacement Estimation . . . . . . . . . . . . . . . . . . . . . 69 5.3.2 Transfer Function Estimation . . . . . . . . . . . . . . . . . . 71 5.3.3 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . 71 5.4 Lateral Relaxation Imaging Results . . . . . . . . . . . . . . . . . . . 74 5.4.1 Displacement Estimation . . . . . . . . . . . . . . . . . . . . . 74 5.4.2 Shear Wave Velocity and Parameter Estimation . . . . . . . . 76 5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.5.1 Axial Relaxation Imaging . . . . . . . . . . . . . . . . . . . . . 77 5.5.2 Lateral Relaxation Imaging . . . . . . . . . . . . . . . . . . . . 79 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.1 Summary and Contributions . . . . . . . . . . . . . . . . . . . . . . . 82 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Appendices A Derivation of an Expression for the Acoustic Radiation Force . . . . . . 97 B The Relationship between the Step and Relaxation Responses . . . . . 101 v List of Tables 2.1 Programmable sequence parameters . . . . . . . . . . . . . . . . . . 20 3.1 Pulse parameters for acoustic output experiment . . . . . . . . . . . 51 3.2 Acoustic output measurements . . . . . . . . . . . . . . . . . . . . . . 53 4.1 Parameter values for ARF experiment . . . . . . . . . . . . . . . . . . 60 5.1 Shear wave velocity and Young’s modulus estimates . . . . . . . . . 77 vi List of Figures 1.1 Ultrasonic imaging system . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 Diagram of main ARF imaging techniques in literature . . . . . . . . 9 2.1 Overall system used in ARF imaging . . . . . . . . . . . . . . . . . . . 16 2.2 ARF Imaging software screen shot . . . . . . . . . . . . . . . . . . . . 18 2.3 ARF pulse sequence with imaging and pushing pulse parameters . 19 2.4 Image formation parameters . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Default TGC curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Comparison of axial and lateral relaxation imaging . . . . . . . . . . 25 2.7 RF echo signal from data import block . . . . . . . . . . . . . . . . . 27 2.8 RF echo signals showing axial displacement . . . . . . . . . . . . . . 29 2.9 Transfer function estimation in axial relaxation imaging . . . . . . . 31 2.10 1D Voigt model used in axial relaxation imaging . . . . . . . . . . . . 33 2.11 Nyquist plot of normalised Voigt model . . . . . . . . . . . . . . . . . 34 2.12 Gelatine phantom used in this work . . . . . . . . . . . . . . . . . . . 42 2.13 Displacement profile showing displacement due to heating . . . . . 44 2.14 Displacement profile showing displacement due to ARF . . . . . . . 45 3.1 Pressure and intensity waveforms with acoustic parameters . . . . . 49 3.2 Axial-lateral scan of pulse intensity integral . . . . . . . . . . . . . . 52 3.3 ARF pulse sequence with imaging and pushing pulse parameters . 56 4.1 Experimental setup for ARF experiments . . . . . . . . . . . . . . . . 59 4.2 Spatiotemporal displacement profile with default parameters . . . . 61 4.3 Results for various numbers of cycles per pushing pulse . . . . . . . 62 4.4 Results for various pushing PRFs . . . . . . . . . . . . . . . . . . . . . 62 vii List of Figures 4.5 Results for various numbers of pushing pulses per sequence loop . 63 4.6 Results for various centre frequencies . . . . . . . . . . . . . . . . . . 63 4.7 Results for various pushing transmit apertures . . . . . . . . . . . . . 64 4.8 Displacement profiles for various pushing transmit apertures . . . . 65 5.1 Spatiotemporal displacement profiles for homogeneous phantoms 70 5.2 Relaxation responses for homogeneous phantoms . . . . . . . . . . 71 5.3 Nyquist plots of relative transfer function estimates . . . . . . . . . . 72 5.4 Magnitude and phase plots of relative transfer function estimates . 72 5.5 Relative elastic and viscous parameter estimates . . . . . . . . . . . 73 5.6 Relaxation time constant estimates . . . . . . . . . . . . . . . . . . . 74 5.7 Lateral relaxation responses for homogeneous phantoms . . . . . . 75 5.8 Time of peak displacement estimation results . . . . . . . . . . . . . 76 5.9 Observed and time constant estimate relaxation responses . . . . . 80 viii Glossary ADC analogue to digital converter ARF acoustic radiation force ARFI acoustic radiation force impulse FDA Food and Drug Administration FIR finite impulse response PRF pulse repetition frequency PRP pulse repetition period PVC polyvinyl chloride RF radio frequency SDK software development kit SDUV shearwave dispersion ultrasound vibrometry SSI supersonic shear imaging SWEI shear wave elasticity imaging TGC time-gain control ix Acknowledgements Firstly, I would like to thank my supervisor, Dr. Tim Salcudean, for his wisdom, guidance and support. He kept me going when nothing worked, and the suc- cesses that resulted owe a lot to his encouragement. The trio of elastography experts, Dr. Ali Baghani, Dr. Hani Eskandari and Dr. Reza Zahiri-Azar, supported me throughout my Masters with excellent ideas, technical help and good humour. The work would have floundered without their assistance. I am very grateful to them all. I would also like to thank the staff at Ultrasonix, particularly Chris Cheung for his help with power supply testing (and letting us melt capacitors off their PCBs), and Corina Leung for her help with acoustic output measurements. Thanks to all the past and present cohabitants of the Robotics and Control Lab for making it a great place to work. The shared ideas, jokes and strange smells will be cherished. Finally, a big thank you to all my friends and the whānau back home, espe- cially Mum, Dad and Haley for always being there. Special thanks to Denver, for making life magic. x —For Nana, suspended in the sky. xi Chapter 1 Introduction 1.1 Motivation and Overview An abundance of ultrasonic imaging techniques exist and their application spans numerous areas. Robotic guidance systems [77], preservation of art works [34], and a host of biomedical procedures make use of ultrasound. Within the biomedical area, ultrasonic imaging is most commonly known for providing im- ages of human foetuses in prenatal care. However, it is used in myriad other ways, such as Doppler imaging of blood flow, high-frequency imaging of the eye [54], and guidance of needle biopsy for prostate and breast cancer. Compared to other common medical imaging modalities, it is portable, inexpensive, safe and fast. An introduction to diagnostic ultrasound imaging and its 60 year history can be read in the text by Szabo [64]. The topic of this thesis is a particular ultrasonic imaging technique named “elastography”, a method of imaging elastic properties of soft tissue. Stiffness is a well-known measure of tissue health, being the parameter “measured” during manual palpation of tissues when inspecting for tumours or other pathological changes in tissue. Palpation, however, does not provide a quantitative measure of stiffness, depends on the person performing it, and is limited to tissues rela- tively close to the surface of the skin. Typical ultrasound B-mode images, based on the amplitudes of ultrasonic echoes returning from the tissue, show contrast between regions of tissue with different acoustic impedance or ultrasonic scat- tering properties, but do not necessarily show contrast between regions with different elastic properties [59]. They therefore sometimes do not provide any detectable contrast between healthy and diseased tissue, particularly in certain cases, such as tumours and liver cirrhosis, where the disease grows out from the healthy tissue [64]. The aim of elastography is to fill this gap by providing 1 1.1. Motivation and Overview images that show the stiffness contrast between healthy and diseased regions of tissue, ideally with quantitative measures that allow for comparison between tissue types and measurement techniques. Such images aid in tissue character- isation and show great potential for assisting with the diagnosis and treatment procedures of cancers of the breast and prostate, for example [65, 83]. These are the most common cancers diagnosed in women and men in Canada in 2010 [6], and it has been shown that such malignant cancers are several times stiffer than healthy tissues [25]. Besides the application to diagnosis, elastography can identify stiffer regions that make suitable targets for biopsy, enhancing the ef- fectiveness of the procedure. Ultrasound elastography has been an active area of research for around 20 years ([46] first uses the term). References [17, 47, 53] provide a good review of elastography techniques and their development. Ultrasonic elastography imaging techniques typically involve scanning a re- gion of the tissue with ultrasound before and sometime after movement has been generated in the medium. Elastic, or viscoelastic, properties of the tissue can be extracted by comparing the state of the medium at the different times. Such techniques require a method of generating movement in the medium, a method of monitoring and measuring the movement, and a method of con- verting the movement into information that can aid in tissue characterisation. Techniques are typically grouped by the method of generating movement. Orig- inally, external, mechanical excitation sources were used [46]; these continue to be in common usage. Natural biological movements created by sources such as cardiac and respiratory cycles also provide a source of excitation which can be used to perform elastography [24]. Another excitation technique that has become a very active research area in the past ten years is the acoustic radiation force (ARF) of ultrasound [40]. Here, nonlinear properties of acoustic wave prop- agation can result in a net force in the direction of the ultrasound propagation [27]. By delivering focused ultrasound to particular regions in tissue, internal, localised movement can be generated. This technique is broadly termed ARF imaging, and many methods exist for implementing and utilising it to perform elastography. ARF imaging is the focus of this thesis. Several methods have been used for monitoring and measuring the tissue movement after excitation, including Doppler-based velocity measurement (e.g. 2 1.1. Motivation and Overview [48, 79]) and displacement or strain estimation using local changes between re- ceived ultrasound echoes with time (e.g. [12, 46]). The latter method is used in the work of this thesis. Non-ultrasound based monitoring has also been used, with systems employing laser vibrometry [69, 70] and magnetic resonance imag- ing [26] both offering alternatives. The final step in elastography, conversion of the monitored movement into useful information regarding elastic properties of the tissue, largely depends on the previous two steps. Raw static strain maps can be used, as stiffer regions of tissue will experience less strain that softer regions, or strains can be con- verted into parameters such as the shear or Young’s modulus, or Poisson’s ratio [47]. If a transient or periodic movement is excited and monitored in the tissue, viscous and frequency-dependent parameters can also be obtained [1, 14, 79]. The impetus behind further exploration of parameters is the potential they pro- vide for improving tissue characterisation. For example, cancerous breast tis- sue can be either benign or malignant. Both types may appear as stiff regions in elasticity-based elastography, but using additional information about viscous and frequency-dependent properties has been shown to aid in separating them, thereby improving tissue characterisation [62]. If an imaging method can pro- vide these additional parameters, the likelihood of correctly diagnosing the type of tumour increases. Returning to ARF imaging in particular, it provides several advantages over external-excitation methods. Firstly, because ARF is generated by a focused ul- trasound transmission, the resulting movement occurs in a relatively localised region (roughly equivalent to the dimensions of the focal region of the trans- mission). This localised nature of ARF imaging means that boundary condi- tions can often be ignored, and that select regions can be palpated indepen- dently. Secondly, the displacements generated by ARF under typical conditions are of the order of tens of micrometres or less. The ultrasonic echoes on which correlation-based displacement estimation is performed are of a high frequency, corresponding to spatial variations of around 100 µm. The low relative displace- ment caused by ARF means that a small search window can be used, such that an incorrect neighbouring peak will not be mistaken for the correct peak (known as “peak hopping”). This results in improved displacement estimation. The small 3 1.2. Background displacement also reduces the stresses imparted on tissue [40]. Finally, the ultra- sound probe used to generate the ARF can also be used to monitor the resulting movement, eliminating the need for external excitation apparatus, and leaving a simplified, integrated system. This thesis presents a new ARF imaging method providing the ability to ex- tract elastic, viscous and frequency-dependent parameters from soft ultrasonic media. The remainder of this chapter contains general background on ultra- sonic imaging and ARF, and a literature review of other ARF imaging techniques, particularly those looking at viscoelastic and frequency-dependent parameters. Chapter 2 presents a full description of the ARF imaging system used in this work, detailing the ultrasound system, associated sequence and processing soft- ware, and the construction of the gelatine tissue phantoms used as the ultra- sonic transmission media in the subsequent experiments. In Chapter 3, issues related to the safety of the methods to the potential patient and to the ultra- sound equipment are studied. Chapter 4 describes initial ARF experiments per- formed to validate the theory presented in the background and test the limits of the ultrasound platform and probe that form the basis of the ARF imaging system. Chapter 5 presents results from experiments with homogeneous tis- sue phantoms performed to test and validate the ARF imaging methods of this work. Finally, the thesis is summarised, conclusion are drawn and suggestions are made for future developments. 1.2 Background 1.2.1 Ultrasonic Imaging This section gives a brief and basic introduction to the principles of pulse-echo ultrasound imaging and phased arrays. For more details refer to the text by Christensen [10]. Ultrasonic imaging takes many forms. However, the most commonly used technique in diagnostic medical imaging applications is pulse-echo ultrasound. Fig. 1.1 presents a typical ultrasonic imaging system that could be used to perform diagnostic medical imaging. The system consists firstly of an ultra- 4 1.2. Background sound platform, which houses the electronic hardware and software required to control the system. This is connected to an ultrasound probe, which con- verts voltage signals produced by the ultrasound platform into ultrasonic pres- sure signals, typically using piezoelectric transducers. In pulse-echo ultrasound, short-duration (order of microseconds), sinusoidal pressure pulses are gener- ated. The surface of the probe is placed in contact with the ultrasonic trans- mission medium, so that, with the aid of acoustic matching layers, the pres- sure pulse conducts into the medium. As the pressure wave of the pulse prop- agates in the medium, echoes are produced at the boundaries between regions with different acoustic impedances. These echoes may travel back through the medium to the probe. Additionally, the cellular makeup of the medium may cause ultrasonic scattering, also producing echo signals that can be received by the probe. The magnitude of the echoes will depend on the acoustic impedances and scattering properties of the different regions of the medium, and the time at which they are received will depend on the distance of these regions from the ultrasound probe. The ultrasound probe converts the echoes back into volt- age signals, which are transferred to the ultrasound platform for processing and display. Essentially, the instantaneous amplitude of the received signal deter- mines the brightness of the pixel corresponding to that time. In this way, a 1D map of depth versus brightness can be displayed, with the brightness giving in- formation about the bulk acoustic and scattering properties with depth. A 2D map, known as a B-mode image (an example of which is shown on the display in Fig. 1.1), can be created by shifting the lateral location at which the transmit- ted pulse is delivered. A common way to shift this lateral location is by using a phased-array ultrasound probe, which contains numerous individual piezoelec- tric transducer elements, any or several of which can be electronically selected by the ultrasound platform. By using multiple elements in each pulse transmis- sion and controlling the relative delay of the pulse applied to the elements, axial and lateral beam focusing and steering can be performed. Fig. 1.1 defines the axial, lateral and elevational directions, which are relative to the probe. These definitions are used throughout this work. In modern ultrasound platforms, the received echo signals are sampled and quantised with an analogue to digital converter (ADC) operating at sampling 5 1.2. Background Figure 1.1: A typical ultrasonic imaging system with the basic components: the ultrasound platform, ultrasound probe, transmission medium and display. The definitions of the axial, lateral and elevational directions are also shown. Ultrasound Platform Transmission Medium Ultrasound Probe Axial Lateral Elevational Display Pulse  6 1.2. Background rates around 20–40 MHz. The echo signals themselves have frequencies around 2–14 MHz, depending on the centre frequency of the elements of the ultrasound probe, and are referred to as radio frequency (RF) echoes because this frequency range is similar to that of electromagnetic RF signals. The digitised RF echoes, or RF data, are typically subject to a number of processing steps, including filtering and envelope detection, before being converted to pixel intensities for display as a B-mode image. Some ultrasound platforms, such as the SonixMDP (Ultrasonix Medical Corp., Richmond, BC) used in this work, have special research modes that allow access to the raw RF data. This allows researchers to explore new tech- niques by which RF data can give useful information regarding the transmission medium. 1.2.2 Fundamentals of ARF The phenomenon of ARF was investigated by Lord Rayleigh in the early 20th century [55, 56], but it was not till the end of the century that it was first used in relation to ultrasound-based tissue characterisation [15, 39, 59, 63]. Even then, it is only at the beginning of the 21st century that this became a very active area of research. The fundamental theory of ARF is well established, but rather complex. Essentially, it is a result of the nonlinear nature of the acoustic wave equation. Considering nonlinear terms results in an “excess pressure” with a non-zero temporal average. This is in “excess” of the usual oscillatory pressure, which has a temporal average of zero. When the wave meets a reflecting or ab- sorbing object, the excess pressure results in a net force, the ARF. Further read- ing is available in [27]. This section outlines some general results applicable to ultrasound-based tissue characterisation. A commonly encountered equation quantifying ARF in soft tissue-like me- dia is F = 2α〈I 〉 c , (1.1) where F is the ARF as a force per unit volume, α is the absorption coefficient of the medium, 〈I 〉 is the temporal average acoustic intensity of the ultrasound and c is the speed of sound in the medium [4, 40, 45]. Appendix A gives a deriva- tion of (1.1). This equation can be applied to each point in the medium, and 7 1.3. Review of ARF Imaging Techniques the force is in the same direction as the ultrasound propagation (i.e. axially). Thus, for a homogeneous medium, the spatial force profile can be considered proportional to the spatial intensity profile. The dependence on 〈I 〉 shows that increases in transmitted pulse amplitude, duration or repetition frequency, all result in increased ARF. The equation assumes plane-wave propagation and that the majority of the attenuation of ultrasound is due to absorption rather than scattering. There are some further nonlinear effects not accounted for by (1.1). At high pressure levels, the density of the transmission medium cannot be considered constant, as the compressional wave creates local variation in density as it prop- agates. Because higher density results in higher wave speed, positive compres- sional pressures travel faster than negative rarefactional pressures. The result is a skewing of the pressure waveform, where positive pressure peaks “catch up” with negative peaks. This effect is used advantageously in harmonic imaging [78], as it creates harmonic frequencies of the fundamental pressure frequency. In ARF imaging, the harmonic frequencies are more readily absorbed by the medium, owing to the frequency dependence of the absorption coefficient [10], and hence the ARF is enhanced. The effect also results in a narrower transmit beam and a slight shift in the acoustic intensity field in the positive axial direc- tion away from the ultrasound probe [40, 58]. 1.3 Review of ARF Imaging Techniques ARF imaging techniques, like all elastography techniques, can be grouped by the method by which excitation is applied to the medium. Typically, three groupings are used: static, transient and periodic. This section examines works from liter- ature covering each of these groups, but pays closer attention to works in which viscoelastic and frequency-dependent parameters are extracted. It also looks at recent trends in ARF imaging, and describes the position of the techniques of this thesis in relation to previous works. A diagram of the works covered and their groupings is shown in Fig. 1.2. 8 1.3. Review of ARF Imaging Techniques Figure 1.2: A diagram showing the groupings by excitation method of important ARF imaging techniques from literature, with references to seminal works for each.  ARF Imaging Static Transient Periodic Acoustic Radiation Force Impulse (ARFI)* McAleavy et al., 2003 Shear Wave Elasticity Imaging (SWEI)* Sarvazyan et al., 1998 Supersonic Shear Imaging (SSI)* Bercoff et al., 2004 Sonorheometry Walker et al., 1997 Acoustic Radiation Force Impulse (ARFI) Nightingale et al . 2001 Harmonic Motion Imaging (HMI) Konofagou & Hynynen, 2003 Vibro-acoustography Fatemi & Greenleaf , 1998 High Frequency Harmonic Imaging Liu & Ebbini 2007 Shear Wave Dispersion Vibrometry (SDUV)* Chen et al., 2009 * indicates method uses monitoring of shear waves 9 1.3. Review of ARF Imaging Techniques 1.3.1 Static Methods In static excitation methods, a fixed stress is applied to the medium and the re- sulting displacement or strain is measured after the system has reached a sta- ble state. Here the frequency of excitation is at or near 0 Hz. Although this is commonly employed in external-excitation elastography, the one main ARF elastography technique using it is acoustic radiation force impulse (ARFI) imag- ing, developed by researchers at Duke University [40, 41, 42]. Here, several long- duration (~20–60 µs), focused pushing pulses are delivered into the tissue, cre- ating an “impulse” in ARF. In the static version of ARFI, the displacement re- sponse at the pushing location is monitored using standard pulse-echo ultra- sound (from the same probe) during and after the pushing pulses, giving an ap- proximate impulse response. Although this would make it a transient method, typical ARFI results are formed only from the displacement at or near the time of peak displacement. In an elastic medium, the peak displacement is inversely related to the Young’s modulus. Therefore, comparison of peak displacements gives a method of presenting relative elastic properties of different media or re- gions of a medium. Other parameters based on the transient response, such as time to peak displacement and recovery time from peak displacement, have also been examined [51]. The same group has developed finite element method simulation techniques for displacement responses and ARF heating effects [49, 50]. ARFI has been applied to ex vivo and in vivo tissues, most recently to the prostate, where it has shown an ability to differentiate between different re- gions of the prostate anatomy, and detect changes due to prostate cancer, be- nign prostate hyperplasia and prostatic calcifications [81, 83]. 1.3.2 Transient Methods Transient ARF methods use the temporal response of displacement or strain to extract properties of the medium in which they are generated. This allows elastic and viscous parameters to be obtained. Transient signals contain a range of fre- quencies, meaning, unlike static methods, transient methods can also be used to examine the frequency dependence of the parameters. One important class of transient (and periodic) methods consists of those 10 1.3. Review of ARF Imaging Techniques that aim to generate and monitor shear wave propagation in the medium. As in ARFI, impulsive pushing pulses are delivered along a certain axial line; how- ever, instead of monitoring the response on this line, the response is measured at several lateral locations some distance from it. The pushing pulses generate an approximately cylindrical shear wave that propagates laterally away from the pushing location, so that the monitored responses show the progression of the shear wave through the medium. Parameters of the shear wave are extracted from these responses and related to viscoelastic parameters of the medium. Shear wave velocity is the most commonly-extracted parameter, as it can eas- ily be related to the shear or Young’s modulus. The idea was first developed as shear wave elasticity imaging (SWEI) by a group based around the University of Michigan and Moscow State University [59]. In SWEI, shear wave velocity is estimated in a homogeneous medium from the change in peak displacement times in responses at different lateral posi- tions. A faster shear wave corresponds to higher elasticity. Later work by related researchers looked at the possibility of extracting viscous parameters from the attenuation of the shear wave as it propagates [18]. The ARFI excitation method has also been used to instigate shear wave prop- agation, with research by the same Duke University group undertaken to de- velop time-to-peak algorithms that were used to obtain shear wave velocity and, from that, shear and Young’s moduli [36, 43, 44]. They have applied their meth- ods in vivo to quantify liver stiffness (with application to detecting fibrosis) [52], looked at effects of viscosity [76] and combined shear wave ARFI with static ARFI [82]. Finally, a group from Université Paris VII developed the supersonic shear imaging (SSI) technique [2], whereby the pushing location is transferred to sev- eral points along the pushing axis. The transition between pushing locations occurs faster than the propagation of the shear wave generated by each push (hence the “supersonic”), and results in a more intense shear wave with a longer axial extent. Using fast ultrasound imaging techniques, a 2D movie of the shear wave propagation is recorded, allowing inversion methods to be used to recon- struct shear elasticity. Since its inception, SSI has been used to quantify breast lesion viscoelasticity in vivo [65], and has been incorporated into a commercial 11 1.3. Review of ARF Imaging Techniques ultrasound imaging platform (SuperSonic Imagine, Aix-en-Provence, France). The same group has also studied viscoelastic and frequency-dependent models of soft tissue [3, 7, 62] and used SSI to extract the parameters of such models [13]. Moving away from transient shear wave methods, early work with transient ARF was performed by a group from the University of Virginia [75]. Eventually termed sonorheometry [74], the method involves transmitting a series of short- duration (< 2 µs) pulses into the medium at a high pulse repetition frequency (PRF). The short duration of the pulses allows their echoes to be processed as imaging pulses, while the high PRF at which they are delivered results in suffi- cient ARF to displace very soft media, such as blood [74] and the vitreous of the eye [73]. The displacement-time response created by the pulses, essentially a step response to the applied ARF, is fitted in the time-domain to a viscoelastic Voigt model, allowing relative elastic and viscous parameters to be obtained, and different media to be compared and classified. A similar technique, though with longer pushing pulses and interspersed imaging pulses, was used by another group [35] in combination with the SWEI imaging technique to give absolute, rather than relative, Voigt model parameters of gelatine tissue phantoms and ex vivo pig muscle. A final transient paper of relevance [48] looks at the relaxation response after removal of ARF of a metal sphere embedded in a homogeneous gelatine tissue phantom. The response is fitted in the time-domain to the theo- retical response of a Voigt model to extract viscoelastic parameters. 1.3.3 Periodic Methods Periodic methods deliver a fixed- or multi-frequency periodic ARF into the medium. The dynamic nature of these methods allows viscoelastic parame- ters to be obtained. If multiple frequencies are used, frequency-dependencies can also be obtained. Like transient ARF techniques, methods may monitor dis- placements at the pushing location or at a lateral location away from the push- ing location to observe the periodic shear waves generated. Vibro-acoustography, by researchers at the Mayo Clinic in Minnesota, was one of the earliest periodic methods developed [8, 15]. Two confocal transducers 12 1.3. Review of ARF Imaging Techniques with transmit centre frequencies separated by a small amount (< 1% of the cen- tre frequencies) interfere to produce an ARF that varies sinusoidally at the differ- ence frequency (an ultrasound equivalent to the beat frequency phenomenon in audio). The difference frequency can be altered, and the amplitude and phase of the periodic response at each frequency can be measured, thereby giving a com- plete frequency response over the swept range. This can then be related to theo- retical responses of viscoelastic models. The group has also experimented with alternative methods of producing a spectrum of excitation frequencies [69, 70] and looked at the shear waves produced. This developed into a technique called shearwave dispersion ultrasound vibrometry (SDUV) [9, 68]. SDUV looks at the dispersion, or frequency dependence, of the shear wave velocity by measuring the phase of the periodic shear wave at different excitation frequencies. Results are fitted to a Voigt dispersion model to extract shear elasticity and viscosity pa- rameters. A similar periodic method is harmonic motion imaging [23], and re- cent work from this group has focused on extraction of viscoelastic parameters and their frequency dependence [72]. Recently, a group from the University of Minnesota began looking at harmonic ARF imaging using high frequency ultra- sound, with application to ophthalmology and thin tissue constructs [28, 29]. Viscoelastic and frequency-dependent parameters are again obtained [19]. 1.3.4 Perspective This review of ARF imaging techniques reveals certain recent trends. Several methods are reaching a mature stage where clinical application is feasible, or realised. An increasing effort is being made to develop a better understanding of the viscoelastic properties of soft tissues, pointed to by the increase in number of works dealing with viscous, not just elastic parameters, and with frequency- dependent changes in these parameters. This latter trend is followed by the work of this thesis. The ARF imaging techniques presented in this thesis, named axial relax- ation imaging and lateral relaxation imaging, use transient excitation methods. In axial relaxation imaging, in contrast to previous excitation methods, push- ing pulses are of medium duration (~7 µs). They are delivered at a high PRFs, 13 1.3. Review of ARF Imaging Techniques creating a step response in the medium. However, instead of monitoring this step response, which is difficult because pushing pulses interfere with imag- ing pulses, the relaxation response after removal of ARF is monitored. The re- laxation response, being a transient signal, has a frequency spectrum contain- ing information about the medium over a range of frequencies. It is analysed in the frequency domain, applying transfer function techniques, and the fre- quency response is related to that of a Voigt model with a frequency-dependent viscous term, governed by a power law. This gives several parameters describing the frequency-dependent viscoelastic properties of the medium, all from a sin- gle relaxation response (acquired in less than 30 ms). Other methods that give this frequency-dependent information typically rely on swept-frequency peri- odic ARF with longer excitation and monitoring times. Lateral relaxation imag- ing, in its currently implementation, is used to provide some validation of the axial relaxation imaging method, and does not differ significantly from previous ARF imaging methods that use shear wave displacement monitored at lateral locations away from the pushing location. 14 Chapter 2 Methodology This chapter outlines the general methods used in this work. It describes the ultrasound equipment, the associated software and the tissue phantoms used in the experiments of later chapters. Two ARF methods were developed in this work, namely axial relaxation imaging and lateral relaxation imaging. Particu- lars of both methods are described in this chapter. An overall diagram of the ARF imaging system used in this work is shown in Fig. 2.1. The left half of the figure consists of the Texo sequence generating software, ultrasound platform and probe, and the transmission medium. The initial generation of ARF and movement, and the imaging of the movement, oc- cur in these blocks. Subsequent processing of data into meaningful information regarding the viscoelastic properties of the transmission medium and the dis- play of this information to the user occurs in the computing device and display blocks in the right half of the figure. 2.1 Ultrasound Platform and Probe The SonixMDP ultrasound platform was used in this work along with the L14–5/38 linear ultrasound probe (Ultrasonix Medical Corp., Richmond, BC). This section outlines the important characteristics of this equipment as it re- lates to ARF imaging. The SonixMDP is an integrated diagnostic ultrasound system with an open software architecture. Through its “Texo” software development kit (SDK) [67], researchers are able to create custom software to control pulse generation and reception. It also gives access to the raw RF echo signals, which can be cap- tured, digitised and saved for later processing. These aspects of the platform are crucial for ARF imaging, as they allow long-duration pushing pulses to be gen- 15 2.1. Ultrasound Platform and Probe Figure 2.1: The overall system used in ARF imaging in this work. Texo Sequence Generating Software Ultrasound Platform Transmission Medium Ultrasound Probe Computing Device Display Data Import Displacement Estimation Axial Relaxation Imaging: Transfer Function Estimation Parameter Estimation Display Conversion Lateral Relaxation Imaging: Shear Wave Velocity Estimation 16 2.2. Texo Sequence Generating Software erated and accurate displacement estimation to be performed on the raw RF echoes. The hardware and software of the SonixMDP delivers the desired elec- trical pulses to the connected ultrasound probe and receives and processes the returned RF echoes. The processing includes filtering, sampling and quantisa- tion of the RF echoes into a form suitable for further digital processing. Control of the SonixMDP through the Texo SDK is discussed in Section 2.2. The L14–5/38 ultrasound probe is composed of 128 elements, arranged in a 1D linear fashion. The elements have a bandwidth from 5–14 MHz with a centre frequency of 7.5 MHz. The ultrasound probe connects directly to the SonixMDP, with 128 lines providing individual connections for each element. The hard- ware of the SonixMDP limits the number of elements that are able to receive RF echoes simultaneously to 64, though it can transmit on up to 128 elements simultaneously. 2.2 Texo Sequence Generating Software The Texo sequence generating software allows the user to control the ultrasound platform. This section describes the software written to perform ARF imaging on the SonixMDP. It also outlines the differences between the sequences used in the axial and lateral relaxation imaging methods. 2.2.1 Texo SDK The sequence generating software is written in the C++ language and employs the Texo SDK [67]. This SDK consists of various C++ classes containing functions and parameters that can be called and modified by the software. These func- tions and parameters control the generation of ultrasound pulse sequences and processing of the received RF echo signals performed by the SonixMDP. Parame- ters such as pulse duration, PRF, transmission voltage and centre frequency can be controlled, as can image formation parameters such as transmit and receive aperture, focal depth, acquisition depth and time-gain control (TGC). The SDK also has a function that allows the digital RF data acquired by the SonixMDP to be saved for later processing. 17 2.2. Texo Sequence Generating Software Figure 2.2: The home screen of the ARF Imaging sequence generating software. 2.2.2 ARF Imaging Software The sequence generating software written for this work, named “ARF Imaging”, provides a command-line interface to the user for controlling the ARF imaging procedures used. It is a heavily modified version of the “console” program pro- vided by Ultrasonix Medical Corp. as a demonstration of the capabilities of Texo. The software only runs on the SonixMDP platform, as it needs to connect to the electronic hardware of the platform in its initialisation. A screen shot of the ARF Imaging home screen is shown in Fig. 2.2. From the home screen, the user firstly needs to select the probe socket (A, B or C) on the SonixMDP to which the desired ultrasound probe is connected. If successful, the program responds with the type of probe connected. The user can then select the type of ARF imaging sequence they want to program, ei- ther axial relaxation imaging or lateral relaxation imaging. After selecting, the program prompts them to either enter values for the sequence parameters or keep the default settings. Most parameters can be modified. The parameters can be broken into three sets: pushing pulse parameters, imaging pulse param- eters and image formation parameters. These are discussed further in Section 2.2.3. After setting the desired sequence parameters, the user can run the se- quence, which begins the transmission of pulses and the sampling of the re- sulting RF echoes. The sequence will loop until terminated by the user (using 18 2.2. Texo Sequence Generating Software Figure 2.3: An example ARF pulse sequence demonstrating the imaging and pushing pulse parameters. .  .  . .  .  . Pushing phaseImaging phase Sequence loop i.1 i.2 i.N i p.1 p.2 p.N p i.1 Imaging pulse Pushing pulse M i  cycles (t PDi ) M p  cycles (t PDp ) V TX T c=1 /f c T PRp=1 /f PRpT PRi=1 /f PRi .  .  . the “Stop sequence” command) or the callback function designed to prevent ex- cessive heating and damage to the probe (see Section 3.3 for more details). The user can then save the acquired RF echo data for later processing by the data processing software. The data are saved in a binary format, with a header that can be read by the data import block of the data processing software. A text file containing the ARF Imaging program version, the sequence type, the ultrasound probe model and values of all the sequence parameters used in the acquisition is also generated. 2.2.3 Sequence Parameters A list of all the sequence parameters, their default values and the ranges over which they can be modified in the ARF Imaging software is given in Table 2.1. Images demonstrating the parameters are shown in Fig. 2.3 and 2.4. Referring firstly to Fig. 2.3, the first part of the sequence shows the imag- ing phase, with the imaging pulse parameters, while the second part shows the pushing phase, with the pushing pulse parameters. The top part of the fig- ure shows individual pulse characteristics, while the lower part shows sequence characteristics. The number of cycles in a pulse, Mi or Mp , is the number of si- 19 2.2. Texo Sequence Generating Software Table 2.1: A list of the sequence parameters that can be programmed with the ARF Imaging software, their default values and the ranges over which they can be modified. Parameter Symbol Default Value Modifiable Range Pushing Pulse No. cycles Mp 48 1–48 PRF fPRp 20 kHz 1–200 kHz Transmit voltage VT X 47 V N/A Centre frequency fc 6.7 MHz 4–20 MHz No. pulses per sequence loop Np 100 1–1000 Imaging Pulse No. cycles Mi 1 1–48 PRF fPRi 10 kHz 1–20 kHz Transmit voltage VT X 47 V N/A Centre frequency fc 6.7 MHz 4–20 MHz No. pulses per sequence loop Ni 300 1–1000 Image Formation Focal depth D f 20 mm 10–50 mm Acquisition depth Da 40 mm 10–50 mm Pushing transmit aperture AT Xp 32 elements 0–128 Imaging transmit aperture AT X i 32 elements 0–128 Imaging receive aperture ARX i 32 elements 0–64 Pushing centre element Xp 64 0–128 Imaging centre element Xi 64 0–128 TGC top gain Gt 40% 0–100% TGC bottom gain Gb 60% 0–100% TGC bottom depth zb D f 10–50 mm TGC middle gain Gm 50% 0–100% TGC middle depth zm 0.5D f 10–50 mm Digital gain GD 2500 0–4095 Receive sampling frequency fs 20 MHz N/A 20 2.2. Texo Sequence Generating Software Figure 2.4: The ultrasound probe and elements, and the transmission medium, shown with the important image formation parameters: focal depth, D f ; ac- quisition depth, Da ; transmit/receive aperture, AT X /RX ; and imaging/pushing centre element, Xi/p . ... ... z A xi al  D ep th Df  ATX/RX Da Lateral Location x 0 0 Xi/p Focus 128 Ultrasound Probe Transmission Medium 1 Probe Elements 21 2.2. Texo Sequence Generating Software nusoidal cycles at the centre frequency, fc , making up the pulse duration, tPDi or tPDp . These parameters are related by tPD = M/ fc , where the appropriate subscript should be applied to tPD and M for imaging or pushing pulses. The transmit voltage, VT X , is the voltage amplitude of the pulse signal generated by the SonixMDP, and is the same for imaging and pushing pulses. The sequence characteristics are described by the PRF, fPRi or fPRp , which is the inverse of the pulse repetition period (PRP), TPRi or TPRp , and by the number of pulses per sequence loop, Ni or Np (labelled in the figure as i.1 to i.Ni for imaging pulses and p.1 to p.Np for pushing pulses). Fig. 2.3 also shows that the sequence loops continuously, with the pushing phase followed by the imaging phase of the next sequence loop. Sequences are set up in this way so that every acquisition begins with an imaging phase, occurring before any pushing pulses are delivered. This provides reference RF echoes for use in displacement estimation (see Section 2.3.2). The pushing phase then creates displacement within the medium, and the subsequent response of the medium is monitored by the imaging phase of the next sequence loop. The sequence effectively creates a relaxation response in the medium: ARF is applied in the pushing phase, then abruptly removed at the beginning of the imaging phase. The ensuing relaxation responses are used to obtain viscoelastic information regarding the medium. Now considering Fig. 2.4, the focal depth, D f , and acquisition depth, Da are shown as axial depths into the transmission medium. The acquisition depth is determined by the time for which RF echo signals are sampled after transmis- sion of an imaging pulse. Time and depth are equivalent here, because of the assumption that the speed of sound in most soft tissue-like transmission me- dia is fixed at 1540 m/s. However, the depth interpretation is more meaningful in this situation, and so is used throughout this work. Typically, the RF echoes from pushing pulses are ignored, as their pulse duration is too long to use in displacement estimation. The parameters AT X i , AT Xp and ARX i (shown collec- tively as AT X /RX in Fig. 2.4) describe the imaging transmit, pushing transmit and imaging receive apertures, respectively. Each is described by the number of elements making up its width. The elements of the L14–5/38 ultrasound probe have a pitch of 0.3048 mm [66], so the number of elements in an aperture can be converted to an equivalent width in millimetres. The centre elements, Xi and Xp 22 2.2. Texo Sequence Generating Software (shown as Xi/p in Fig. 2.4), define the lateral location of the centre of the relevant aperture. The lateral location of the pushing centre element is typically used as the 0 mm reference for the lateral axis. This is mainly of importance in lateral relaxation imaging, where the lateral locations of imaging centre elements are set at some distance from the pushing location. Several parameters are not shown in Fig. 2.3 and 2.4, as these are related to the pre-processing performed by the SonixMDP on the received RF echo sig- nals. The digital gain, GD , simply amplifies the RF echoes to a suitable level before they are sampled and quantised by the ADC in the SonixMDP. The value can be programmed between 0 and 4095. The approximate mid-point, 2500, was found to give good results in most experiments. The receive sampling fre- quency, fs , is the sampling frequency at which the RF echoes are sampled by the ADC. The TGC curve consists of five parameters that control the way in which the RF echo signals are amplified as a function of depth (or, equivalently, time). This goal of this function is to account for the attenuation of the medium, so that deeper echoes, which have been subject to more attenuation, are amplified more. These five TGC parameters are the gain at zero depth (the “top” of the curve) , Gt ; the gain and position at the maximum depth (the “bottom” of the curve), Gb and zb , respectively; and the gain and depth of a mid-point between the top and bottom,Gm and zm , respectively. The gain values are numbers from 0–100%, representing percentages of the full gain range of the integrated circuit that performs the TGC function in the SonixMDP. The default value of zb is set to the value of the focal depth, D f , while the default value of zm is half the value of zb . Note that depths beyond the maximum depth, zb , have the same gain, Gb , applied to them. Together, the parameters describe a piecewise linear TGC curve, the default version of which is shown in Fig. 2.5. 2.2.4 Axial and Lateral Relaxation Imaging Sequences As mentioned, two sequence types have been created in this work to perform ARF imaging. These are termed axial relaxation imaging and lateral relaxation imaging. Both types employ the “push-then-image” technique described in Sec- tion 2.2.3. In the pushing phase of this technique, higher-intensity pushing 23 2.2. Texo Sequence Generating Software Figure 2.5: The default TGC curve showing the five parameters: top gain, Gt , bottom gain, Gb , bottom depth, zb , middle gain, Gm , and middle depth, zm . 0 100 50 G ai n (% ) Depth (mm) 40302010 zb zm Gm Gb Gt pulses are focused at a specific axial depth and lateral location to create local movement of the medium through ARF. The difference between axial and lat- eral relaxation imaging comes in the imaging phase. In axial relaxation imaging, the lower-intensity imaging pulses monitor the axial movement at the same lat- eral location as the pushing pulses (i.e. on the same axis). In lateral relaxation imaging, the imaging pulses monitor the axial movement at a lateral location some distance from the original pushing location. Fig. 2.6 shows the difference in imaging location between axial and lateral relaxation imaging. In axial relaxation imaging, a single acquisition is sufficient to provide vis- coelastic information about the medium at the position at which the acquisition is performed. This is because the relaxation response at this position is alone used to perform force-displacement transfer function estimation, as outlined in Section 2.3.3. In lateral relaxation imaging, two or more acquisitions from different lateral locations are required to perform shear wave velocity estimation, as outlined in Section 2.3.5. This is typically done by running the acquisitions separately. Each acquisition is performed with a different lateral imaging location, entered manually by the user. Here, all other sequence parameters are fixed and the ul- trasound probe is not moved between acquisitions. Another method, which can 24 2.2. Texo Sequence Generating Software Figure 2.6: A comparison of axial and lateral relaxation imaging, showing the different imaging locations used in each sequence type. Ultrasound Probe Transmission Medium Axial Relaxation Imaging x z Pushing Lateral Relaxation Imaging Push creates movement Movement imaged at same location Movement imaged at different location Pushing location Imaging location be done with the ARF Imaging software, is to image multiple lateral locations in alternating sequence loops of one acquisition. Here, a single lateral location is imaged after a given pushing phase, and a different lateral location is imaged af- ter the next pushing phase. This method, however, can cause issues if the move- ments from a previous push have not completely died out, as they will affect the response monitored at the next lateral location. A final method is to alternate the lateral locations with each imaging pulse, so that each sequence loop monitors the movement from a single push at a variety of locations. The disadvantage of this approach is that the imaging PRF needs to be increased by a factor equal to the number of imaging locations, which may strain the limits of the ultrasound system. The first method (multiple, separate acquisitions) is used in this work for its trouble-free simplicity. To move the lateral imaging location, the imaging centre element, Xi , is adjusted for each acquisition. Each output spatiotempo- ral profile therefore shows how the axial displacements change with time at a fixed lateral location. The lateral imaging locations are spaced by approximately 1 mm or more, so there is an appreciable difference between the resulting pro- files. Each location is also approximately 2 mm or more from the lateral location at which the pushing pulses are delivered, as results are unpredictable within the region of excitation of the pushing pulses [44]. The lateral dimension of the re- 25 2.3. Data Processing Software gion of excitation can be approximated by the lateral beam width, which, for the L14–5/38 ultrasound probe and the default pushing pulse parameters of Table 2.1, is approximately 0.94 mm [10]. 2.3 Data Processing Software The data processing software consists entirely of MATLAB (The MathWorks, Natick, MA) m-code running on a PC. MATLAB is a mathematics analysis pack- age suitable for large-scale data processing, and its m-code allows the user to create routines and functions that execute MATLAB functions. The digitised RF echo data obtained by the SonixMDP are transferred offline to a PC. The pro- cessing of these RF echoes can be broken into several blocks, as shown in Fig. 2.1. These blocks are termed: data import, displacement estimation, transfer function estimation (for axial relaxation imaging), shear wave velocity estima- tion (for lateral relaxation imaging), parameter estimation and display conver- sion. All blocks are implemented in MATLAB m-code. This section describes each of the blocks. 2.3.1 Data Import The digitised RF echo data from the imaging device are imported into MATLAB using the data import block1. This involves reading the header of the data file, which specifies information regarding the acquired data: the number of samples in each RF echo, the data format of each sample, the total number of RF echoes, the number of RF echoes per sequence loop, and the total number of sequence loops. Optionally, the number of sequence loops imported can be limited, to reduce processing time. The block then organises the data after the header into a matrix with each column corresponding to an individual echo signal. Further processing can then be performed on this matrix. A typical RF echo signal output by the data import block is shown in Fig. 2.7. The full output from the block typically consists of several thousand of these 1The m-code for this block is based on m-code by Sergei Koptenko (Resonant Medical, Mon- treal, QC) written in February 2005. 26 2.3. Data Processing Software Figure 2.7: A typical RF echo signal as output by the data import block. 0 5 10 15 20 25 30 35 40 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 104 RF echo amplitude with depth Axial Depth (mm) A m p l i t u d e  echoes. The amplitude is given as a quantisation count, and does not have an associated unit. 2.3.2 Displacement Estimation The displacement estimation block uses the changes between RF echo signals to estimate the local 1D axial displacement along the acquisition depth2. It uses time-domain correlation with sub-sample cosine interpolation [12]. The displacement estimation algorithm uses a reference RF echo signal from a particular lateral location. This is typically acquired when the transmission medium is at rest (i.e. before applying any pushing pulses), thus giving a snap- shot of the 1D speckle pattern of the medium prior to motion. Subsequent RF echoes from the same lateral location are compared with this reference to esti- mate the motion that has occurred since the reference signal was acquired. All displacement estimates are therefore relative to the reference signal. Because the echoes are acquired at a regular rate, the motion can be recorded as a func- 2The m-code for this block is based on m-code by Reza Zahiri-Azar (University of British Columbia, Vancouver, BC) written in October 2008. 27 2.3. Data Processing Software tion of time, with a sampling rate equal to the imaging PRF, fPRi . The individual RF echo signals can be up-sampled from the receive sample frequency, fs , be- fore displacement estimation. This improves the accuracy of the displacement estimates. In the processing of this work, all echoes were up-sampled by a factor of ten. A reasonably high factor was used because data were processed offline and processing time was not critical. To estimate the axial motion at multiple depths, each RF echo is divided into a number sections, or windows. The size of each window determines the degree of locality of the displacement estimates, and can be set within the m-code. A value of 1 mm was used in the processing performed in this work. Each window from the reference RF echo is compared with the equivalent window from the RF echo for which displacements are being estimated (the “current” window). The reference window is incrementally shifted across the current window in a nor- malised 1D cross-correlation. The range over which it is shifted is defined by the search radius, which should be set greater than the peak displacement expected. It can be set in the m-code, and a value of 100 µm was used in the processing of this work. The point over this range at which the correlation score is maximum gives the displacement estimate, with sub-sample cosine interpolation used to increase the accuracy of this estimate. Fig. 2.8 shows an example of equivalent parts of two RF echo signals taken from the same location in a transmission medium at different times. It can be seen in certain features (marked A, B and C), that the current echo pattern has displaced in relation to the reference echo. The displacements of the features are marked as uA , uB and uC . Although the displacement estimation block does not rely on matching of point features, by considering the features as centres of windows, the figure serves to illustrate the concept that different axial depths can move by different amounts. The m-code outputs a spatiotemporal displacement profile as a matrix, giv- ing the axial displacement as a function of depth into the medium and time. An example of such a profile can be seen in Fig. 2.14, later in this chapter. 28 2.3. Data Processing Software Figure 2.8: An example of two RF echo signals from the same location in a trans- mission medium at different times. Current RF echo Reference RF echo A A B B C C uA uB uC z Axial Depth (mm) 18.5 19 19.5 20 20.5 29 2.3. Data Processing Software 2.3.3 Axial Relaxation Imaging: Transfer Function Estimation This section describes the transfer function estimation that follows displace- ment estimation of RF echo signals acquired with axial relaxation imaging. The processing is performed on the spatiotemporal displacement profile output by the displacement estimation block. The transfer function estimated for axial relaxation imaging is that between the force applied at the focal depth (i.e. the ARF) and the resulting displacement at the focal depth. The relationship between force and displacement (or, equiv- alently, stress and strain) is fundamental to the viscoelastic behaviour of the medium. Therefore, the transfer function of this relationship is closely linked to the viscoelastic properties. The method assumes that the system relating the force, f (t ), and axial dis- placement, u(t ), at the focal depth (z = D f ) is linear, and that the force can be treated as a scaled, negative step function f (t )=K (1−θ (t )) , (2.1) where θ (t ) is the Heaviside step function and K is an unknown, constant scal- ing factor that depends on the magnitude of the ARF. This assumption implies the force reaches a constant level over each set of pushing pulses, which, as dis- cussed in Section 4.4, is valid for the typical parameters used in this work. It also implies the force drops to zero immediately following cessation of the push- ing phase, which seems a reasonable assumption given that the pushing pulses stop and the imaging pulses that begin are of very low intensity in compari- son. Although K is unknown, the focal-depth displacement relaxation response, u (t )= r (t ), to f (t ) for t ≥ 0, can be converted to a system step response relative to K and used to obtain a relative transfer function. The relative step response, sr (t ), is calculated as sr (t )=K s (t )= r (0)− r (t ) , (2.2) where s (t ) is the unit step response to f (t ) = θ (t ). This expression is valid for causal, linear, time-invariant systems that do to respond instantaneously (i.e. the impulse response is zero for t ≤ 0). See Appendix B for a proof of this. Be- 30 2.3. Data Processing Software Figure 2.9: The system diagram for identifying the force-displacement transfer function in axial relaxation imaging. Transfer Function H(ω) Displacement u(t) Force f(t) cause the derivative of a step function is an impulse function, the assumed lin- earity of the system allows the relative impulse response, hr (t ), to be obtained from the relative step response by differentiation in time. Thus, hr (t )=Kh (t )= d dt sr (t )=− d dt r (t ) , (2.3) because r (0) is a constant. Here, h (t ) is the (non-scaled) impulse response. Fi- nally, taking the Fourier transform of (2.3) gives the relative, force-displacement transfer function estimate, Hr (ω)=KH (ω)=K U (ω) F (ω) =F {hr (t )} , (2.4) where ω is angular frequency,U (ω) and F (ω) are the Fourier transforms of u (t ) and f (t ), respectively, andF {·} represents the Fourier transform operation. The diagram of Fig. 2.9 shows the system diagram used for identifying the force- displacement transfer function in axial relaxation imaging, including the as- sumed form of the force, f (t ). The MATLAB m-code of this transfer function estimation block first obtains the relaxation response at the focal depth, r (t ), by taking the mean of the re- sponses over a small axial range centred at the focal depth. This helps reduce noise in the relaxation response while still preserving the primary information at the focal depth. The relative impulse response, hr (t ), is obtained accord- ing to (2.3) using the diff derivative function in MATLAB. The relative transfer function is then calculated with the fft fast Fourier transform function. The re- 31 2.3. Data Processing Software sulting transfer function estimate consists of a vector of complex numbers, each of which corresponds to the value of Hr (ω) at a number of discrete frequencies. The frequencies depend on the imaging PRF, fPRi , and the number of sampled time instances in the displacement relaxation response, which is equal to the number of imaging pulses per sequence loop, Ni . The transfer function esti- mate is input to the parameter estimation block, and processed as described in the next section. 2.3.4 Axial Relaxation Imaging: Parameter Estimation Estimates of the parameters of a viscoelastic model of the transmission medium can be obtained by relating the model to the transfer function estimate. This processing does not necessarily need to be employed, as simple parameters of the transfer function itself can be used directly without parameter estimation. If parameter estimation is used, numerous methods are possible, as several differ- ent viscoelastic models are available, and there are various ways to relate these to the estimated transfer functions. However, in this work, the transfer functions estimated through axial relaxation imaging are related to the 1D Voigt model. The 1D Voigt model consists of a linear elastic element, µ, in parallel with a lin- ear viscous element, η, as shown in Fig. 2.10. The relationship between stress, σ, and strain, ε, for each of these elements is σe =µε (2.5) and σv = ηε̇, (2.6) where the subscripts e and v denote the elastic and viscous components, re- spectively. The complete stress-strain relationship for the model is then σ=σe +σv =µε+ηε̇. (2.7) Taking the Fourier transform of (2.7), the stress-strain transfer function for the 32 2.3. Data Processing Software Figure 2.10: The 1D Voigt model used in axial relaxation imaging, shown in rela- tion to the transmission medium. µη Probe Transmission Medium Voigt Model Focal Point 33 2.3. Data Processing Software Figure 2.11: Nyquist plot of normalised Voigt model transfer function, HV (ω). 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 Real I m a g i n a r y Nyquist plot of normalised Voigt model  1D Voigt model is HV (ω)= ε̂ (ω) σ̂ (ω) = 1 µ+ jωη , (2.8) where ε̂ (ω) and σ̂ (ω) are the Fourier transforms of ε and σ, respectively. To give an idea of the form of this model, the Nyquist plot of HV (ω) for values ofµ= 1 Pa and η= 1 Pa·s is shown in Fig. 2.11. In this plot, frequency increases as the curve is followed from right to left, beginning withω= 0 at the right-most point (where HV (ω)= 1). The curve describes a semi-circle, and is characteristic of first-order linear systems such as the 1D Voigt model [11]. It can also be shown that the magnitude response of the Voigt model has a cut-off frequency, fcutoff = 1 2pi µ η , (2.9) and that the phase response is ∠HV (ω)= tan−1 (−ωη µ ) . (2.10) Considering the proportional relationship in 1D between stress and force, 34 2.3. Data Processing Software and strain and displacement, the transfer function of (2.8) can be equated with that of (2.4) using a new constant of proportionality, K ′, say. Thus Hr (ω)=K ′ 1 µ+ jωη , (2.11) where Hr (ω) is the relative transfer function estimated in axial relaxation imag- ing. Equation (2.11) allows the parameter estimation processing block to obtain relative values of the viscoelastic parameters of the Voigt model. The relative elastic parameter, µr , can be found as the inverse of the DC component of the relative transfer function. I.e. at ω= 0, µr = µ K ′ = 1 Hr (ω= 0) . (2.12) The relative viscous parameter, ηr , is related to the imaginary component of the relative transfer function: ηr = η K ′ = 1 ω ℑ{H∗r (ω)} |Hr (ω)|2 , (2.13) where ℑ {·} denotes the imaginary part of the argument, and the superscript ∗ stands for complex conjugation. The relative viscous parameter can theoretically be evaluated at any value of ω > 0 for which an estimate of Hr (ω) exists. There are two issues with this, however. Firstly, inaccuracies will be present in the estimate of Hr (ω), particu- larly where its frequency components are of very low magnitude and dominated by noise. This effect can be mitigated by using frequency values with higher magnitude. The second issue is due to the limitations of the Voigt model in predicting the dynamic behaviour of the transmission medium. Gelatine and soft tissue media have been shown to exhibit frequency-dependent variations in model parameters than cannot be accounted for by a simple 1D Voigt model. Thus, a different value of ηr will be obtained at each value of ω. Models that ac- count for these effects quickly become complex. However, some authors [14, 22] have noted that a power-law relationship between the viscous parameter and 35 2.3. Data Processing Software frequency may be used with a Voigt model to give a measure of the frequency dependence of the viscous parameter. Such a model treats this relationship as η (ω)= ηcωn , (2.14) whereηc is referred to as the consistency index andn as the flow index. Typically, n < 0 in gelatine and soft tissue media, implying that viscosity decreases with frequency, as is often observed [3, 7, 14]. The relationship of (2.14) can be fitted to the value of ηr versus ω obtained by (2.13) to give an estimate for the flow index. The consistency index can only be estimated to a relative value, ηc,r = ηc/K ′, due to the relative nature of ηr , and thus the flow index, n, is the real parameter of interest in this model. A final parameter that is investigated in this work is the Voigt model time constant, τ, defined as τ (ω)= η (ω) µ = ηr (ω) µr . (2.15) This parameter is of interest because, in theory, it is not dependent on the mag- nitude of the ARF (i.e. it is not a relative parameter) but still conveys information regarding frequency-dependent properties of the medium. In fact, because µ is a constant, the same flow index, n, can be obtained from τ (ω). The MATLAB m-code implementing this parameter estimation method uses the relationships of (2.12), (2.13) and (2.14) to obtain estimates for the parame- ters µr , ηr and n, respectively. The vector of complex numbers that constitutes Hr (ω) and the associated vector of frequencies that forms ω are used directly in equations (2.12) and (2.13), where basic MATLAB functions are employed. A lin- ear least-squares regression of the log of (2.14) is used to recover the flow index, which is obtained as the gradient of the regression line. 2.3.5 Lateral Relaxation Imaging: ShearWave Velocity Estimation Lateral relaxation imaging works by monitoring the shear wave that propagates laterally from the pushing location after removal of ARF at the end of the pushing phase. As a shear wave propagates, it is influenced by the transmission medium. In particular, elastic properties of the medium tend to determine the propaga- 36 2.3. Data Processing Software tion velocity of the shear wave, while viscous properties tend to attenuate it. By monitoring two or more different lateral locations, these changes in the shear wave can be observed. In this work, the shear wave velocity is estimated us- ing the time of peak displacement in the responses at different lateral locations (the lateral relaxation response). This method or similar has been applied in nu- merous previous works (see Section 1.3), and is used in this work primarily as a means to provide some validation of the results obtained with the axial relax- ation imaging method. Shear wave velocity estimation begins with estimation of the time of peak displacement in the focal-depth relaxation response of each lateral imaging lo- cation acquired. These estimates are denoted tpk1, tpk2, . . . , tpkQ , where Q is the total number of lateral relaxation responses acquired and x1,x2, . . . ,xQ are their lateral locations in millimetres. The estimation can be performed manually by inspecting plots of each relaxation response, or automatically, using MATLAB to find the peak value of the response and the corresponding time index. The lat- ter method is used in this work. Each relaxation response is pre-processed as in axial relaxation imaging, by taking the mean of the responses over a small axial range centred at the focal depth. The result is filtered with a tenth-order, low- pass, zero-phase finite impulse response (FIR) filter with a cut-off frequency of 500 Hz to remove jitter on the displacement estimates. Using a zero-phase fil- ter ensures no delay is added to the responses that could affect estimation of the time of peak displacement. The max function in MATLAB is then used to find the index of the peak displacement, so that the value of tpk can be obtained from the element of the time vector with the same index. A plot is then made for each phantom with the values of x along the x-axis and the corresponding values of tpk along the y-axis, and a linear least-squares regression is performed. The inverse of the gradient of this regression line is taken as the shear wave velocity estimate, cs . If responses were acquired in only two locations (Q = 2), then the estimate can be simply calculated as cs = x2−x1 tpk2− tpk1 , (2.16) assuming x1 is closer to the pushing location than x2. 37 2.3. Data Processing Software In estimating cs , some consideration may need to be made as to whether or not the medium is homogeneous. If not, the lateral extent of the values of x used in the regression should be limited to a suitable pixel size, over which the medium can be treated as homogeneous. 2.3.6 Lateral Relaxation Imaging: Parameter Estimation After obtaining a shear wave velocity estimate with lateral relaxation imaging, it can be related to the Young’s modulus, E , of the medium using the relationship cs = √ E 2(1+ν)ρ , (2.17) where ν is the Poisson’s ratio of the medium and ρ is the density [36]. For gela- tine and soft tissue media, it can generally be assumed ν= 0.5 (implying incom- pressibility) and ρ = 1000 kg/m3, giving the equation E = 3000c2s . (2.18) Note that the lateral relaxation imaging method used in this work does not currently give viscous or frequency-dependent parameter information. 2.3.7 Display Conversion The display conversion block uses one or more of the quantities (or derived quantities) from either the transfer function estimation block, the shear wave velocity estimation block or the parameter estimation block, and converts it into a form suitable for the display. It uses suitable MATLAB commands, such as plot and imagesc to send the information to the display. In this work, the con- version and display occur offline after acquisition of the RF echo signals. How- ever, potential exists for processing and display to occur in near real time, pos- sibly in conjunction with other ultrasound imaging modes, such as B-mode. 38 2.4. Transmission Medium 2.4 TransmissionMedium 2.4.1 Background The transmission medium in Fig. 2.1 is intended to represent any material in which ARF-induced displacements can be generated and ultrasonically tracked, such as human or animal soft tissue, or tissue-mimicking phantoms. Phan- toms are made from materials such as gelatine or polyvinyl chloride (PVC). They are designed to replicate the acoustic and mechanical properties of soft tissues, particularly the speed of sound, acoustic absorption and scattering, and elastic and viscous parameters. In this work, experiments are performed with gelatine- based phantoms. These phantoms are described in this section. Gelatine has long been used in the construction of tissue-mimicking phan- toms for ultrasound, because its mass density is close to 1000 kg/m3, its speed of sound closely matches that of human soft tissues, and additives can be used to control the acoustic attenuation [30, 32]. Graphite, glass beads or cellulose are typically used as additives. The individual particles in these substances are tens of micrometres in diameter, and when distributed throughout a gelatine mould, produce an ultrasound speckle pattern similar to that seen in soft tis- sues. Increasing the concentration of the additive increases the acoustic at- tenuation (both absorption and scattering) [30, 57]. This is important for ARF imaging, as the absorption of the medium in part determines the magnitude of the ARF, and the speckle pattern is necessary for displacement estimation. The speed of sound can be controlled through the addition of n-propanol or glyc- erol to achieve a value very close to the commonly accepted soft tissue value of 1540 m/s [30, 31, 32, 57]. The mechanical properties of gelatine have been studied, and it has been shown to be a suitable material for mimicking the viscoelastic behaviour of soft tissues [20, 33]. Increasing the gelatine concentration increases the elasticity and viscosity of the resulting phantom, though at concentrations less than 10%, the elasticity has a greater dependence on gelatine concentration than viscosity [14, 21]. The Young’s modulus has been measured for various gelatine concen- trations and it has been shown by several authors that there is an approximate 39 2.4. Transmission Medium square-law relationship, reported empirically in [20] as Eg = 0.0034C2.09g , (2.19) where Eg is the Young’s modulus (in kPa) of the gelatine phantoms in that work and Cg is the concentration (in g/L). Note that measurement of Eg can be vary with boundary conditions and test methods. Nonetheless, (2.19) is used for pre- dicting elastic values of the phantoms in this work. Gelatine phantoms are commonly employed in ARF works. The Duke Uni- versity group created very soft gelatine phantoms using the method of [20], pro- ducing phantoms with values of Eg from 0.1–1.6 kPa [40]. This method is also used in [28]. The sonorheometry group from the University of Virginia has de- veloped very soft acrylamide-based phantoms to model the acoustic and me- chanical properties of the vitreous of the eye [38]. However, these phantoms are not stiff enough for the intended applications of this work. 2.4.2 PhantomConstruction The phantoms used in this work were made with ~300 bloom, type A, porcine- skin gelatine (Sigma-Aldrich Corp., St Louis, MO) and type 50 Sigmacell cellu- lose (Sigma-Aldrich Corp.). This cellulose has a 50 µm average particle size, and provides absorption and scattering. These dry ingredients were hydrated with distilled water. The following process was used to construct the phantoms: 1. Heat the distilled water to 70°C in a suitable beaker on a hot plate. 2. Add the gelatine and cellulose to the water and stir continuously. 3. Bring the solution to 80°C, then remove heat. Continue stirring continu- ously. 4. Allow the solution to cool to 30°C, then pour it into the desired container. 5. Cover the container and place in a refrigerator until congealed (overnight, typically). 40 2.5. Development Issues A precision balance (BH-3000, Excell Precision Co., Ltd., Santa Clara, CA) was used to measure the correct quantities needed to achieve the desired concen- trations. A stirring hot plate (PC-420D, Corning Inc., Lowell, MA) was used for heating and stirring, along with a magnetic stirring bar (Fisher Scientific, Ot- tawa, ON). The speed of sound in phantoms made with this recipe should be very close to 1540 m/s, though a small variation will likely occur with gelatine and cellulose concentrations [20, 57]. The viscoelastic parameters will depend on the gelatine concentration, while the attenuation (and hence ARF magnitude) will depend on the cellulose concentration. All parameters are dependent on temperature, which should be controlled in experiments. Note that because no anti-fungal agent was added, the life of the phantoms constructed with this process is lim- ited to a few weeks, if they are kept refrigerated when not in use. A small quantity of bleach could be used to extend their life [37], however this was not necessary for the work performed here. The gelatine concentrations of the phantoms produced were between 1.5 and 4 wt% and the cellulose concentration was always 4 wt%. The low gela- tine concentrations produced phantoms that were very soft. This was required to ensure significant displacements could be generated, as the SonixMDP pro- duces a relatively weak ARF push (see Section 2.5 for more information). The low gelatine concentrations also allowed the cellulose to settle somewhat to the bottom of the container. Despite this, the cellulose concentration in the bulk of the phantom provided sufficient attenuation to generate ARF displacements. All phantoms were made with 200 mL of distilled water and set in clear plas- tic cups, as shown in Fig. 2.12. The height was approximately 80 mm and the diameter at the top surface of the phantom was 65 mm. These were covered and kept refrigerated when not in use. 2.5 Development Issues Numerous obstacles were faced in the development of the techniques described in this chapter. This section highlights some of these, and explains ways around them. The aim of including this information is to prevent the same issues oc- 41 2.5. Development Issues Figure 2.12: A gelatine phantom used in this work. curring for others who may want to follow parts of this work. 2.5.1 Heating and ARF Although ultrasound attenuation creates ARF in a transmission medium, it also causes heating. This is discussed in relation to safety issues in Section 3.1.1, but it also has effects on the ARF imaging methods, particularly the displacement estimation step. The two primary effects of heating are thermal expansion and the change in the speed of sound in the medium with temperature. In linear thermal expansion, a material that is heated will expand according to its coefficient of thermal expansion. In soft tissues undergoing ARF imaging, this coefficient is typically so low and the temperature increases so small that effects of thermal expansion can be ignored. Simulations by the Duke Univer- sity group of their ARFI imaging method [49] showed that displacements due to thermal expansion are very small; 0.02 µm was the largest axial displacement magnitude they observed. The change in the speed of sound, also simulated in [49], has a more signifi- cant effect. As the medium heats during ultrasound imaging, the speed of sound may increase or decrease, causing objects to appear nearer or further, respec- tively, from where they were originally. Over time, they seem to displace towards 42 2.5. Development Issues or away from the ultrasound probe. The change can be modelled with a linear coefficient, kc , which can be positive or negative, depending on the medium. I.e. ∆c = kc∆T, (2.20) where ∆c is the change in the speed of sound, kc is the linear coefficient and ∆T is the change in temperature, which should be less than 6°C for the linear assumption to hold [49]. Values of the linear coefficient have been reported for breast fat, kc = −2.93 m/s/°C, cow liver, kc = 1.48 m/s/°C, and a rubber phan- tom, kc =−3.8 m/s/°C [61]. In water and most tissues, kc is positive, though it is negative for fatty tissues [60]. The simulations of [49] gave a maximum displace- ment error generated by an ARFI imaging sequence of approximately 2.5 µm, using the above kc value for liver. They noted, however, than such an error had not been observed experimentally and that their simulations of the temperature increases contained several factors that tended to over-estimate the error. These studies were important in the development of the techniques of this work. Early in the development, when imaging phantoms with short-duration, high-PRF ultrasound and processing the RF echoes with displacement estima- tion, displacements in the medium were observed. These displacements oc- curred over several seconds, were always away from the ultrasound probe and their overall magnitude increased with PRF. It was initially unclear from where the displacements originated, and ARF seemed a potential source. However, as shown in Fig. 2.13, the spatiotemporal displacement profiles were puzzling, as there was no distinct peak at the focal depth, and the displacements did not ap- pear to saturate, even after a minute of imaging. It was not until understanding the concepts of thermal expansion and, particularly, the change in the speed of sound that can occur with heating, that heating was seen as a more likely expla- nation. The fact that all displacements at all depths were directed away from the ultrasound probe indicated the speed of sound in the phantom was decreasing with temperature, or had a negative value of kc . The displacements could not have been generated by ARF, because of the lack of a focal peak and saturation. The time scale was closer to that observed in heating than ARF, where displace- ments occur over milliseconds. These experiments were performed on fairly stiff 43 2.5. Development Issues Figure 2.13: The spatiotemporal displacement profile (left) and the focal depth displacement time series (right) in a firm phantom where displacements are due to heating effects, not ARF. Time (s) isplacement (µm)  D e p t h  ( m m ) D  0 2 4 6 8 10 5 10 15 20 25 30 35 -10 0 10 20 30 40 50 60 60 0 2 4 6 8 10 12 -10 0 10 20 30 40 50 Time (s) D i s p l a c e m e n t  ( µ m ) Displacement time series at depth of 20 mm  PVC phantoms, and no ARF displacements were being generated. This was re- vealed when ARF displacements were finally generated in a very soft, 1.5 wt% gelatine phantom, as shown in Fig. 2.14. Here, a strong peak in displacement at the focal depth and a fast increase in displacement, which saturated after ap- proximately 10 ms, was observed. 2.5.2 Generating Pushing Pulses Another major issue in developing ARF imaging methods concerned the gener- ation of suitably strong pushing pulses. The SonixRP ultrasound platform (Ul- trasonix Medical Corp.) was initially used with a Texo parameter called trex, in an attempt to produce long-duration, continuous-wave, ARF pushing pulses. This parameter allows pulse durations of up to 500 µs [80]. In Texo, trex can be toggled on or off whenever a continuous-wave pushing pulse is required. Unfor- tunately, this switching caused unpredictable corruption of the subsequent RF echo signals used to image the medium. One to four RF signals would be cor- rupted, either consisting only of zero readings, showing large oscillations like the continuous-wave ARF pushing pulses themselves, or showing large ampli- tude spikes throughout the signal. An attempt was made to insert “dummy” imaging pulses between the continuous-wave pushing pulse and the first de- 44 2.5. Development Issues Figure 2.14: A spatiotemporal displacement profile (left) showing ARF- generated displacements in a very soft gelatine phantom. An obvious peak oc- curs at the focal depth of 20 mm, the displacement time series for which is also shown (right). Time (s) isplacement (µm)  D e p t h  ( m m ) D  0 2 4 6 8 x 10-3 0 5 10 15 20 25 30 35 40 -1 0 1 2 3 4 5 5 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 D i s p l a c e m e n t  ( µ m ) Time (s) Displacement time series at depth of 20 mm  sired imaging pulse. These pulses were intended to absorb the corrupted sig- nals until the system corrected itself. However, despite the addition of up to four dummy pulses, corruptions in subsequent pulses still occurred. Because of the issues with the trex parameter, an attempt was made to in- stead generate ARF displacements using a series of high-PRF pushing pulses, rather that a single continuous-wave pulse. Using an oscilloscope to observe the output voltage waveforms, it was discovered that the pulse duration of stan- dard pulses (i.e. those generated with trex disabled) was limited to 1.2 µs. Addi- tionally, changing the PRF or pulse duration during a sequence resulted in cor- ruption of RF echoes similar to that which was seen with the trex parameter enabled. Despite numerous experiment with Texo parameter settings, these is- sues were only resolved when the newer SonixMDP platform was used instead of the SonixRP. Although the trex parameter was still unreliable, corrupted signals were no longer observed when changing pulse parameters within a sequence. It was also discovered that the duration of a standard pulse could be programmed up to 48 cycles, or 7.2 µs at a 6.7 MHz transmit centre frequency, greatly aiding ARF generation. There remains, however, an issue with the power supply architecture used 45 2.5. Development Issues in the Ultrasonix ultrasound platforms. When it experiences high power draw, as occurs during series of high-PRF, long-duration pushing pulses, the voltage level progressively drops from its nominal level of 47 V [80]. The lowered voltage reduces the power in the generated pulses, limiting the magnitude of the ARF and preventing stiffer media from being displaced. The voltage also takes time to recover to its original level after the high power draw stops, and this affects imaging pulses that occur after a pushing pulse. The problem is more severe on the SonixRP than the SonixMDP, but even with the latter, data issues are expe- rienced when exceeding a PRF of approximately 20 kHz or increasing the push- ing transmit aperture beyond 32 elements. The problem was investigated in ex- periments at Ultrasonix Medical Corp., and appears to lie with the rating of the capacitors that provide instantaneous power to the elements of the ultrasound probe [16]. 46 Chapter 3 Safety Issues This chapter examines some of the important safety issues related to ultrasound-based ARF imaging. It relates important acoustic output parameters of ARF imaging to established regulatory limits for standard pulse-echo ultra- sound imaging. An experiment to measure these parameters for the ARF imag- ing method of this work is described. The chapter concludes with a discussion of heating issues in ARF imaging. 3.1 Background Several biological effects occur in human tissue insonified by standard pulse- echo ultrasound imaging [64]. These can be separated into thermal and me- chanical effects, with the most important effects of each being heating and cavi- tation. These effects are addressed by different parameters and limits defined in the U.S. Food and Drug Administration (FDA) guidance document for diagnos- tic ultrasound imaging systems and transducers [71]. 3.1.1 Heating Heating of tissue undergoing ultrasound examination occurs due to conduction and absorption. Conduction is typically from the ultrasound probe surface into the tissue, as the elements of the probe heat during imaging. Absorption oc- curs within the tissue as the ultrasound wave travels through it. Heating due to conduction from the ultrasound probe is typically not an issue, as good probe- tissue coupling means most of the delivered energy is transmitted into the tissue as an ultrasonic wave, and blood perfusion reduces any heating that does occur. The primary heating mechanism in diagnostic ultrasound imaging is therefore 47 3.1. Background absorption. Such heating is of concern, as temperatures in excess of 39°C have the potential to damage tissue [64]. The amount of heating depends on the ab- sorption properties of the tissue and the time-averaged intensity of the deliv- ered ultrasound pulses. The spatial-peak, temporal-average intensity parame- ter, ISPTA , and its associated limit in the FDA guidance document, is designed to prevent harm due to heating from occurring to a patient. 3.1.2 Cavitation Although absorption in tissue results in ARF, a mechanical effect, the resulting movement is typically very small and will not cause harm. Thus, the primary mechanical effect of concern is cavitation. Cavitation occurs when ultrasonic pressure waves create gas bubbles within the medium which subsequently col- lapse or resonate, releasing significant energy and potentially damaging the tis- sue in which they formed. The threshold for the occurrence of cavitation is a function of the peak rarefactional pressure, pr , and centre frequency, fc , of the delivered ultrasound pulse [64]. The mechanical index parameter, MI , a unit- less value calculated from pr and fc , is given a limit in the FDA guidance doc- ument to ensure cavitation does not occur in a patient undergoing diagnostic ultrasound. 3.1.3 Acoustic Output Parameters The previous sections introduced the important parameters, ISPTA and MI , from the FDA guidance document. Before delving further into these parameters, it is worth considering the waveform of acoustic pressure and intensity typically used by the pulse sequences of ARF imaging. An example of such a waveform is shown in Fig. 3.1. It is assumed that these waveforms represent the pressure and intensity at the spatial location at which the various parameters are maxi- mum. Each pulse can be described by a pressure-time waveform consisting of a number of sinusoidal cycles of frequency fc , with a total pulse duration of tPD, as shown in the upper part of Fig. 3.1. The time between each pulse is the PRP, TPR , the inverse of the PRF, fPR . The peak rarefactional pressure, pr , is impor- tant for determining the mechanical index, MI . The intensity waveform, shown 48 3.1. Background Figure 3.1: The pressure and intensity waveforms for a fixed-PRF pulse se- quence, shown with the important acoustic parameters from the FDA guidance document. Pressure Intensity t PD T PR PII I SPPA  I SPTA p r in the lower part of Fig. 3.1, is derived from the pressure waveform by I = p 2 ρc , (3.1) where I is the instantaneous intensity, p is the instantaneous pressure and ρ and c are the density and speed of sound of the transmission medium, respectively [10]. The time-integral of the intensity over the pulse duration, shown by the shaded area in Fig. 3.1, is the pulse intensity integral, PI I . The time-average of the intensity over the pulse duration, shown by the upper dotted line, is the spatial-peak, pulse-average intensity, ISPPA . The time-average of the intensity over the PRP, shown by the lower dotted line, is the spatial-peak, time-average intensity, ISPTA . Values of pr , ISPPA and ISPTA are typically measured in a water tank envi- ronment. However, delivering equivalent ultrasound pulses into tissue will yield lower values, due to the significantly higher attenuation of tissue, which dissi- pates more pulse energy as they travel to the spatial peak. To account for this effect, the FDA guidance document specifies a standard derating that can be ap- plied, which is based on an assumed ultrasonic attenuation coefficient in tissue 49 3.2. Acoustic Output Experiment of 0.3 dB/cm·MHz in the axial direction. The subscript “.3” is added to derated parameters. For example, the derated peak rarefactional pressure is [64] pr.3 = pr exp (−0.0345 fczSP ) , (3.2) where zSP is the axial coordinate of the spatial location of the peak rarefactional pressure. Note that in the exponent, fc must be expressed in MHz and zSP in cm. Similarly, applying the derating to the spatial-peak, time-average intensity gives [64] ISPTA.3 = ISPTA exp (−0.0691 fczSP ) , (3.3) where the exponential attenuation factor doubles because intensity is propor- tional to the square of pressure. The mechanical index can then be calculated as [71] MI = pr.3√ fc . (3.4) The intensity parameters are all interrelated [64]: PI I = ISPPA·tPD = ISPTA·TPR . (3.5) The limits imposed by the FDA are MI < 1.9, ISPPA.3 < 190 W/cm2 and ISPTA.3 < 720 mW/cm2. Note that these limits are the derated values and apply only to “peripheral vessel” ultrasound imaging. Different limits are in place for cardiac, fetal and other applications. 3.2 Acoustic Output Experiment An experiment was performed to determine the acoustic output delivered by the ARF sequences used in this work and compare it with the limits defined by the FDA guidance document [71] and those of other ARF works [2, 40, 42, 65]. 50 3.2. Acoustic Output Experiment Table 3.1: The values of the pulse parameters used in the acoustic output exper- iment. Parameter Value Pulse duration, tPD 0.42 µs PRF, fPR 5 kHz Transmit voltage, VT X 47 V Centre frequency, fc 9.5 MHz Focal depth, D f 20 mm Transmit aperture, AT X 32 elements 3.2.1 Method A water tank acoustic intensity measurement system with a membrane hy- drophone (Onda Corporation, Sunnyvale, CA), situated at Ultrasonix Medical Corp., was used to take measurements of a SonixMDP ultrasound platform with an L14–5/38 linear probe connected. A custom sequence was created in Texo to control the parameters of the delivered pulses. Initially, a very high pulse length and PRF were used. However, due to the need to run the sequence continuously for several minutes to obtain the required measurements, this damaged the el- ements of the probe before the measurements could be made. To avoid this problem, a shorter pulse length and PRF were used, and the results were inter- polated to the longer pulse length and higher PRF typically used in ARF imaging. After replacing the damaged probe, the pulse parameters given in Table 3.1 were set with Texo. These parameters allowed the sequence to be run continuously without any risk of damaging the probe. The interpolation was performed by scaling the measured value of PI I by the ratio of the pushing or imaging pulse duration to the experimental pulse duration (0.42 µs). Equation (3.5) was then used to calculate values of ISPPA and ISPTA for each pulse duration. These values were derated using the appropriate centre frequency for the pulse (6.7 MHz for pushing and imaging pulses and 9.5 MHz for the experiment pulses). In this work, a default ARF pushing pulse is 7.2 µs in duration (48 cycles at 6.7 MHz) and is delivered at a PRF of 20 kHz. This means the pulse duration is 17.1 times that used in the acoustic output 51 3.2. Acoustic Output Experiment Figure 3.2: An axial-lateral 2D scan of pulse intensity integral, PI I , measured in a water tank during the acoustic output experiment. Pulse Intensity Integral Lateral Location (mm) A xi al  D ep th  (m m ) m J/ cm 2 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 65 55 45 35 25 15 70 60 50 40 30 20 10 -8 -4 0 4 8-1 0 -6 -2 2 6 10 experiment, and the experimental value of PI I can be scaled by this factor to give the interpolated pushing pulse value. A typical imaging pulse in this work has a pulse duration of 0.15 µs (1 cycle at 6.7 MHz) and a PRF of 10 kHz, so that the experimental value ofPI I can be scaled by 0.36 to give the interpolated value for a typical imaging pulse. 3.2.2 Results The hydrophone was scanned spatially in the water tank to find the location of peak acoustic output. After finding the peak elevational intensity, a 2D axial- lateral scan was performed at this elevation. Measurements of PI I were taken at each scanning step to give a 2D map, shown in Fig. 3.2. The peak in PI I is centred laterally (0 mm) and occurs at an axial depth of 20 mm, corresponding with the centre element and focal depth of the sequence parameters used. After finding the spatial peak output, the water tank measurement system was used to obtain spatial peak measurements of PI I , MI and ISPPA . These 52 3.2. Acoustic Output Experiment Table 3.2: The acoustic output measurements for the low-power pulses used in the experiment and the interpolated values for ARF pushing and imaging pulses. Pulse Type MI PI I (mJ/cm2) ISPPA.3 (W/cm2) ISPTA.3 (mW/cm2) Experiment 0.41 0.096 61.2 129 Pushing 0.48 1.64 90.1 13000 Imaging 0.48 0.034 90.1 135 values and the interpolated values for typical ARF pushing and imaging pulses are given in Table 3.2. Note that the value of MI was adjusted using (3.4) to account for the different centre frequencies used in the experiment and in the ARF sequences. 3.2.3 Discussion The experimental, pushing and imaging pulse values of MI are well below the limit of 1.9. Similarly, the values of ISPPA.3 are well below the 190 W/cm2 limit. These results imply that cavitation will not occur in tissue undergoing ARF imag- ing with the system of this work, as the instantaneous pressures generated are not high enough. Clearly, the interpolated value of ISPTA.3 for the pushing pulses exceeds the 720 mW/cm2 limit. It should be noted, however, that the true value is likely to be considerably smaller, as the power supply issues discussed in Section 2.5.2 will cause the pressure amplitude to drop over both the duration of the pulse and over multiple pulses. Nonetheless, the very high value, 13000 mW/cm2, is 18 times the limit, and certainly of concern. Although it seems to indicate that ARF imaging will cause dangerous heating in tissue, the relationship between the parameter and the heating it causes needs to be considered further. The limit imposed by the FDA on ISPTA.3 assumes that the pulsing sequence runs indefinitely. This may be the case in standard diagnostic imaging, but in ARF imaging, only a limited number of long-duration pushing pulses are delivered before the regular, short-duration imaging pulses recommence. The pulsing se- quence is thus more complex, and it becomes difficult to apply the FDA limits of [71] to this case, for which the guidance document was not designed. 53 3.2. Acoustic Output Experiment In the literature, works regarding the ARFI imaging technique report values for ISPTA.3 of a similar order to the results here [40, 42]. For example, they mea- sure a sequence with an 29 V transmit voltage, 2.2 µs pulse duration, 7.2 MHz centre frequency and 3.6 kHz PRF to have an ISPTA.3 of 2400 mW/cm2 [40]. They note here that although the value exceeds regulatory limits, the application of this time-average intensity lasts for only a few milliseconds, the time it takes to complete a scan of a region, and therefore is unlikely to generate temperature increases above 1°C, making it safe for human application. Their later simula- tions confirm this is the case for most situations [49]. Works regarding the SSI technique [2, 65], approach the issue of the applica- bility of the ISPTA.3 limit to ARF imaging differently. In the method of [2], four or five pushing pulses are delivered at a PRF of 500 Hz, followed by very low in- tensity imaging pulses. In assessing the safety of the method, they effectively calculate the PI I for a single pushing pulse, then determine how many can be applied per second without violating the 720 mW/cm2 limit. This seems to be a reasonable approach, as it effectively gives a longer time-average over which ISPTA.3 is calculated. Applying this concept to typical sequences in this work, which use 100 push- ing pulses and 300 imaging pulses per sequence loop, the values of PI I in Table 3.2 can be used to show that ten sequence loops per second could be performed while still meeting the limit on ISPTA.3. Although this argument may be reason- able, it must be treated with caution. For example, it could be argued that over ten seconds, 100 sequence loops could be run without violating the FDA limit. However, if those acquisitions were made over the first second, followed by nine seconds without pulsing, there obviously could still be significant risk of over- heating tissue in that first second. There ultimately are no regulatory guidelines tailored to the variable pulsing schemes of ARF imaging. Despite this, it seems reasonable to conclude that so long as the number of acquisitions per second is kept at or below ten, and these acquisitions are spread evenly in time, the tech- niques of this work are safe for use on human tissue. 54 3.3. Probe Heating and Protection 3.3 Probe Heating and Protection As noted in Section 3.2, if an ARF sequence is run continuously, the elements of the ultrasound probe will be permanently damaged. Damaged probes are not only expensive to replace, but the heating that causes the damage could be a safety issue if the probe is in contact with a patient. To prevent such heating, a method was developed to limit the energy delivered from the ultrasound plat- form to the probe. The energy delivered to each element depends on the pulse transmit volt- age, duration and PRF, and the total time for which a sequence is run. As seen in Section 2.2.3, the ARF sequences used in this work are somewhat complex, with different pulse durations and PRFs for pushing and imaging pulses. Fig. 2.3 (reprinted here as Fig. 3.3 for reference) shows the generalised ARF sequence, and is used for the computations in this section. Here, each sequence loop con- sists of an imaging phase followed by a pushing phase. The imaging phase is a series of Ni imaging pulses (labelled i.1 to i.Ni ) transmitted at the imaging PRP, TPRi . The pushing phase is a series of Np pushing pulses (labelled p.1 to p.Np ) transmitted at the pushing PRP, TPRp . Each imaging pulse consists of Mi cy- cles at the centre frequency, fc , while each pushing pulse consists of Mp cycles at the same frequency. Both imaging and pushing pulses have transmit voltage amplitude VT X . In the ARF Imaging Texo program, described in Section 2.2, the sequence loop is repeated until terminated by the user or a command within the program. The energy delivered to the transducer, ET , over L repetitions of the sequence, assuming a 50W probe impedance, is given by ET = V 2T X 100 MiNi +MpNp fc L. (3.6) This quantity represents the power delivered by the transmitted signal (V 2T X /100) multiplied by the total “on” time of pulses. Note that, although not explicit in (3.6), there is a dependence on PRP that arises from the relationship between L and the total time for which the sequence 55 3.3. Probe Heating and Protection Figure 3.3: An example ARF pulse sequence demonstrating the imaging and pushing pulse parameters. .  .  . .  .  . Pushing phaseImaging phase Sequence loop i.1 i.2 i.N i p.1 p.2 p.N p i.1 Imaging pulse Pushing pulse M i  cycles (t PDi ) M p  cycles (t PDp ) V TX T c=1 /f c T PRp=1 /f PRpT PRi=1 /f PRi is looped, tT . That is, tT = ( NiTPRi +NpTPRp ) L. (3.7) Thus, a sequence that employs shorter imaging and pushing PRPs (or higher PRFs) will deliver energy more quickly, implying it is of a higher power. This raises the question of whether an energy or a power limit is required to prevent damage to an ultrasound probe. The heating of a transducer within an ultra- sound probe will begin with the starting of the sequence. If the power of the sequence is low enough, the heating of the probe will stabilise, due to its abil- ity to dissipate heat into the surrounding environment. However, if the power is too high, the heating will cause the temperature to continue to rise to a point where the transducer is damaged. Whether or not this point is reached, depends on the energy delivered to the probe. Because the sequences in this ARF work typically use high PRFs, long pulse durations and high transmit voltages, it can be assumed in most cases that if it were run long enough, the power of the se- quence is high enough to cause damage to the probe. With this assumption, a limit on the total delivered energy seems more appropriate. If it is assumed that a maximum delivered energy, ET,MAX , is allowable be- 56 3.3. Probe Heating and Protection fore damage to the transducer occurs, then the maximum number of sequence repetitions, LMAX , can be determined from (3.6). I.e. LMAX = ET,MAX 100 V 2T X fc MiNi +MpNp . (3.8) A suitable value for ET,MAX for the SonixMDP platform with an L14–5/38 probe was found empirically by considering an early ARF experiment in which a high power sequence was applied to an ultrasound probe for a period of time without causing any damage. In this experiment, the sequence parameters were L = 35,VT X = 47 V, fc = 6.7 MHz, Mi = 1, Ni = 300, Mp = 48, and Np = 100, giving a total energy, ET = 0.589 J. ConsideringVT X will remain at 47 V in all ARF exper- iments of this work, (3.8) can be re-written in a form suitable for programming into the ARF Imaging program: LMAX < fc 38 ( MiNi +MpNp ) , (3.9) where fc is in Hz. In the Texo program, the value of LMAX is calculated when the sequence is first created. When acquisition is started by the user, a function is called after each sequence loop is acquired. The function increments a counter and checks whether LMAX has been reached. If it reaches this number before the user stops acquisition, the program automatically stops acquisition, preventing damage to the probe. The user can also manually reduce the number of sequence loops with the Texo program. This is recommended if only a few sequence loops are required, as it will ensure minimal heating of the probe. The user can then save the acquired data. The empirical sequence repetition limit is known to prevent damage, but it is not known how much higher the limit could be set. If required, experiments could be performed to determine this. However, in this work, the limit is suit- able, and serves to prevent accidental damage to the probe. Note that the theory here could be applied to other ultrasound equipment, but a different value for ET,MAX would need to be determined. 57 Chapter 4 ARF Experiments 4.1 Introduction This chapter details experiments performed to confirm some of the theory re- lated to ARF imaging, verify that the ultrasound system was capable of generat- ing ARF-induced displacements in soft tissue-like transmission media, and ex- plore the limits of these capabilities. It examines the effects of certain pushing pulse and image formation parameters. All experiments used the axial relaxation imaging sequence described in Sec- tion 2.2.3 and a gelatine phantom constructed with the method described in Section 2.4.2. Displacement relaxation responses were generated and imaged in the phantoms, and the peak values of these responses were used to give a mea- sure of the magnitude of the ARF created by the sequence parameters. 4.2 Method The experiments of this section were performed with the SonixMDP platform and an L14–5/38 ultrasound probe. The transmission medium was a gelatine phantom made with 3 wt% gelatine and 4 wt%, type 50 cellulose. The phantom was left unrefrigerated overnight to reach room temperature (21°C). Distilled water was poured into the phantom cup so that a layer approximately 10 mm deep sat on the top surface of the phantom. The ultrasound probe was clamped to a metal strut on an anti-vibration table, and its face was dipped just below the water line, so that pulses were transmitted through approximately 7 mm of water to the phantom. This setup is shown in Fig. 4.1. In each experiment, the parameter of interest was varied using the ARF Imaging Texo software, while all other parameters were kept the same at their 58 4.2. Method Figure 4.1: The experimental setup used for the experiments of this section. The ultrasound probe is dipped into the water layer on top of the phantom, and clamped to a metal strut on an anti-vibration table. default values. The default parameter values given in Table 2.1 were used for all sequence parameters, except the one being varied in a given experiment. Note that this means the lateral pushing and imaging locations were the same, so the sequences were the same as used in axial relaxation imaging. The param- eters examined were: the number of cycles per pushing pulse, Mp ; the pushing PRF, fPRp ; the number of pushing pulses per sequence loop, Np ; the centre fre- quency, fc ; and the pushing transmit aperture, AT Xp . Table 4.1 outlines the val- ues that were tested for each parameter. Note that the values selected for fc are the only valid values over this range, due to the frequency-division method by which the SonixMDP generates the transmit centre frequency [67]. To minimise variations due to atmospheric conditions, all acquisitions in the experiment for a given parameter were taken within eight minutes or less, and all acquisitions over all parameters were taken on the same day over approximately three hours. The acquired RF echoes were processed with the data import and displace- 59 4.3. Results Table 4.1: The values used for each parameter in the ARF experiment. Parameter Values Mp 12, 24, 36, 48 cycles fPRp 10, 20, 30,..., 80 kHz Np 25, 50, 75, 100, 125, 150 pulses fc 4.0, 5.0, 6.7, 8.0, 10.0, 13.3 MHz AT Xp 16, 32, 48, 64 elements ment estimation blocks of the data processing software. A window size of 1 mm and search radius of 100 µm were used in the displacement estimation for all acquisitions. The generated spatiotemporal displacement profiles were used to give relaxation responses at the 20 mm focal depth. These responses were formed by averaging over a depth of 2 mm (2 windows) at the focal depth. For a given parameter, the relaxation responses for each tested value of that param- eter were plotted together, and the peak values of the responses (at t = 0) were plotted against the parameter values. 4.3 Results The spatiotemporal displacement profile for the default parameter values (Mp = 48, fPRp = 20 MHz, Np = 100, fc = 6.7 MHz and AT Xp = 32) is shown in Fig. 4.2. The profile shows that the peak displacement occurs at the beginning of the imaging phase, t = 0, and at the focal depth, D f = 20 mm, and that the overall displacements exhibit a relaxation response, decreasing in magnitude with time. The results for each tested parameter – Mp , fPRp , Np , fc and AT Xp – are shown in Figs. 4.3–4.7. In all results, excepting those where significant displace- ment was not observed (e.g. the fc results for 4 and 13.3 MHz), it was confirmed that the displacement peak indeed occurred at the focal depth. In these results also, the correlation coefficients of focal-depth displacement estimates were all greater than 0.9. Note that to better show the individual relaxation responses on the left of each figure, they are shown only over the first 5 ms, even though the imaging phase lasted for 30 ms. Also, the result for fPRp = 70 kHz was omitted, 60 4.4. Discussion Figure 4.2: The spatiotemporal displacement profile (left) generated in a 3 wt% gelatine phantom using default sequence parameters, and the relaxation re- sponse at the focal depth (right). Time (ms) lacement (µm)  D e p t h  ( m m ) Disp  0 1 2 3 4 10 15 20 25 30 35 5 10 15 20 25 30 0 2 4 6 8 10 12 14 16 18 20 0 5 10 15 20 25 Time (ms) D i s p l a c e m e n t  ( µ m ) Relaxation response at focal depth  as it was very inconsistent with the results for other values, showing very little displacement in comparison. A final figure, Fig. 4.8, shows the spatiotemporal displacement profiles for the various values of AT Xp tested. A change in the size of the region in which displacements are produced can be seen with the changing aperture. 4.4 Discussion The spatiotemporal displacement profile obtained using default sequence pa- rameters, on the left in Fig. 4.2, shows the basic form of the relaxation response obtained in all results of this chapter. The peak displacement, approximately 28 µm, occurs at the focal depth, 20 mm, and at the beginning of the imaging phase, t = 0. The subsequent displacement relaxation appears to follow an ex- ponential curve, decreasing towards 0 µm, as seen clearly on the right in Fig. 4.2. The nature of this response depends on the viscoelastic properties of the medium in which it occurs, and this principle is employed in axial relaxation imaging. The results in Fig. 4.3, for Mp , the number of cycles per pushing pulse, show that peak displacement, and therefore ARF, increases with Mp . The trend is ini- 61 4.4. Discussion Figure 4.3: The results for the experiment in which the number of cycles per pushing pulse, Mp is varied. The focal depth relaxation response for each value tested is shown on the left, while the peak displacement versus the tested value is shown on the right. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 5 10 15 20 25 30 Time (ms) D i s p l a c e m e n t  ( µ m ) Relaxation responses at focal depth for various values of M  p  12 cycl 30 es 24 cycles 36 cycles 48 cycles 10 15 20 25 30 35 40 45 50 5 10 15 20 25 Number of cycles per pushing pulse, Mp P e a k d i s p l a c e m e n t  ( µ m ) Peak displacements for various values of Mp  Figure 4.4: The results for the experiment in which the pushing PRF, fPRp , is varied. The focal depth relaxation response for each value tested is shown on the left, while the peak displacement versus the tested value is shown on the right. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 5 10 15 20 25 30 35 40 Time (ms) D i s p l a c e m e n t  ( µ m ) Relaxation responses at focal depth for various values of pushin  g PRF  38 10 20 30 40 50 60 70 80 18 20 22 24 26 28 30 32 34 36 Pushing PRF (kHz) P e a k d i s p l a c e m e n t  ( µ m ) Peak displacements for various values of pushing PRF 10 kHz 20 kHz 30 kHz 40 kHz 50 kHz 60 kHz 80 kHz  62 4.4. Discussion Figure 4.5: The results for the experiment in which the number of pushing pulses per sequence loop, Np , is varied. The focal depth relaxation response for each value tested is shown on the left, while the peak displacement versus the tested value is shown on the right. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 5 10 15 20 25 30 Time (ms) D i s p l a c e m e n t  ( µ m ) Relaxation responses at focal depth for various values of N  p  25 pul 29 ses 50 pulses 75 pulses 100 pulses 125 pulses 150 pulses 20 40 60 80 100 120 140 160 20 21 22 23 24 25 26 27 28 Pushing pulses per sequence loop, Np P e a k d i s p l a c e m e n t  ( µ m ) Peak displacements for various values of Np  Figure 4.6: The results for the experiment in which the transmit centre fre- quency, fc , is varied. The focal depth relaxation response for each value tested is shown on the left, while the peak displacement versus the tested value is shown on the right. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -5 0 5 10 15 20 25 30 35 Time (ms) D i s p l a c e m e n t  ( µ m ) Relaxation responses at focal depth for various values of centre frequen  cy  4 MHz 5 MHz 6.7 MHz 8 MHz 10 MHz 13.3 MHz 4 5 6 7 8 9 10 11 12 13 0 5 10 15 20 25 30 Centre frequency, fc (MHz) P e a k d i s p l a c e m e n t  ( µ m ) Peak displacements for various values of centre frequency  63 4.4. Discussion Figure 4.7: The results for the experiment in which the pushing transmit aper- ture, AT Xp , is varied. The focal depth relaxation response for each value tested is shown on the left, while the peak displacement versus the tested value is shown on the right. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 5 10 15 20 25 30 Time (ms) D i s p l a c e m e n t  ( µ m ) Relaxation responses at focal depth for various values of pushing trans  mit aperture  16 elem 30 ents 32 elements 48 elements 64 elements 15 20 25 30 35 40 45 50 55 60 65 12 14 16 18 20 22 24 26 28 Pushing transmit aperture (no. elements) P e a k d i s p l a c e m e n t  ( µ m ) Peak displacements for various values of pushing transmit aperture  tially linear, with the displacement doubling from 9 to 18 µm as Mp doubles from 12 to 24. This agrees with the basic theory of ARF, particularly (1.1). The time- average acoustic intensity in this equation, 〈I 〉, is proportional to Mp (Mp essen- tially controls the duty cycle of the pushing pulses), and the ARF magnitude, F , is proportional to 〈I 〉. Assuming the phantom can be treated as a viscoelastic Voigt medium and the displacement reaches a steady state over the pushing phase in each experiment, the steady-state displacement, assumed to be equal to the peak displacement in the results of this chapter, is proportional to F . Hence, the peak displacement should be proportional to Mp . After these initial results, the peak displacement increases more slowly as Mp increases. This is likely to be due to saturation in the power that can be generated by the SonixMDP, with the power supply issues discussed in Section 2.5.2 resulting in reduction of the transmit voltage with large Mp . This effect seems to be apparent in the other results of this chapter also. In any case, the maximum value, Mp = 48 cycles, which is the maximum allowed by the SonixMDP firmware, gives the largest dis- placement, and therefore is the most suitable for ARF generation. In a similar result is observed in Fig. 4.4 for pushing PRF, fPRp . Peak dis- placement generally increases with fPRp , but the rate of that increase reduces 64 4.4. Discussion Figure 4.8: The spatiotemporal displacement profiles for each tested value of pushing transmit aperture, AT Xp . Time (ms) D e p t h  ( m m ) Displacement (µm) for Ap = 16 elements   0 5 10 15 10 15 20 25 30 35 0 5 10 15 Time (ms) D e p t h  ( m m ) Displacement (µm) for Ap = 48 elements   0 5 10 15 10 15 20 25 30 35 0 5 10 15 20 Time (ms) D e p t h  ( m m ) Displacement (µm) for Ap =   32 elements  0 5 10 10 15 20 25 30 35 15 0 5 10 15 20 25 Time (ms) D e p t h  ( m m ) Displacement (µm) for Ap =   64 elements  0 5 10 10 15 20 25 30 35 15 0 2 4 6 8 10 12 14  65 4.4. Discussion at higher values of fPRp . The peak displacement eventually saturates at 60 kHz, and even reduces at 80 kHz. Following the argument for Mp above, time-average intensity, ARF and peak displacement should all be proportional to fPRp , be- cause fPRp also controls the pushing pulse duty cycle. Again, the high levels of power drawn from the SonixMDP power supply at high values of fPRp likely cause the transmit voltage to drop and limit the peak displacement that can be generated. The 20 kHz default value typically used in the experiments of this work generates sufficient ARF in most cases, but it can be seen from these results that an increase by a factor of around 1.4 could be achieved if it were increased to 60 kHz. This is not a huge increase, given the additional heating the higher PRF will generate in the ultrasound probe, but may be useful in some cases. Unlike Mp and fPRp , the number of pushing pulses per sequence loop, Np , does not affect the time-average intensity and the magnitude of the ARF gen- erated. Instead, it affects the total time for which the ARF is applied. This al- lows the effect of displacement saturation due to the elasticity of the medium to be investigated. If a constant force is applied to a viscoelastic medium, the displacement will follow an exponential curve, initially increasing quickly, then slowing its rate of increase till it reaches an asymptote as it saturates. The plot of peak displacement versus Np in Fig. 4.5 shows that the peak displacement gradually saturates over a range from 25–100 pulses per sequence loop, match- ing well the displacement saturation that occurs over time in an elastic medium experiencing a constant force. This was also observed in the right-hand plot of Fig. 2.14. This means that though larger values of Np increase the peak displace- ment, eventually increasing Np will have little effect. The reduction in peak dis- placement seen after Np = 100 in Fig. 4.5 is again likely due to the power supply issues in the SonixMDP. Thus, Np = 100 is an optimal default value, as it en- sures saturation of displacement without over-straining the power supply. Dis- placement saturation is an assumption made in the transfer function estimation method of axial relaxation imaging (Section 2.3.3). The results for the centre frequency, fc , experiments essentially reveal the bandwidth of the L14–5/38 ultrasound probe. This probe has a nominal band- width from 5–14 MHz, with a centre frequency at 7.5 MHz, though information from Ultrasonix Medical Corp. indicates the upper bandwidth limit is somewhat 66 4.4. Discussion lower than 14 MHz [66]. Fig. 4.6 shows that no noticeable displacements were generated at centre frequencies of 4 and 13.3 MHz, outside or near the upper limit of the bandwidth of the probe. Between these frequencies, the plot follows a bandwidth-type curve, peaking at 8 MHz, right near the centre frequency of the probe. The slower roll-off in peak displacement at frequencies above 8 MHz may be related to the nonlinear effect whereby ARF is enhanced at higher fre- quencies, due to the higher absorption of the medium (see Section 1.2.2). The results clearly show that 6.7 or 8 MHz should be used with this probe to ensure best ARF generation. Finally, the peak-displacement results for the pushing transmit aperture, AT Xp , shown in Fig. 4.7, show that an aperture of 32 elements gives the most displacement. It gives a 28 µm peak displacement, twice that generated by 16 elements. It is again likely that the power supply issues discussed in Section 2.5.2 reduce the ARF and peak displacement at higher apertures, due to the ad- ditional power draw of the added elements. Thus, AT Xp = 32 elements is an optimal value. The spatiotemporal displacement profiles of Fig. 4.8 show that as the aperture widens, the axial extent of the region over which displacements are observed increases. This concurs with the concept of f-number, defined as the ratio of the focal length to the aperture width. A wider aperture reduces the f-number, which narrows the depth of focus, or, in the case of focused ul- trasound, the axial extent of the region of excitation. Based on the 0.3048 mm element pitch in the L14–5/38 probe, the f-numbers for apertures of 16, 32, 48 and 64 with a 20 mm focal depth are 4.1, 2.1, 1.4 and 1.0, respectively. The results of this chapter confirm the basic theory regarding ARF imaging, as presented in Section 1.2.2. They also provide a lot of useful information re- garding the capabilities and limitations of the SonixMDP ultrasound platform, and guide the choice of sequence parameters used in subsequent work. 67 Chapter 5 Homogeneous Phantom Experiments 5.1 Introduction The axial and lateral relaxation imaging techniques were applied to homoge- neous phantoms to asses their ability to provide viscoelastic information about the phantoms. Data were acquired for several such phantoms, and transfer function estimation, shear wave velocity estimation, and parameter estimation were performed for each. To provide some validation, the axial relaxation imag- ing results were compared with empirical predictions and results given by the lateral relaxation imaging technique. This chapter first outlines the details of the method used to perform the experiment. It gives the results of axial relaxation imaging then lateral relaxation imaging, and, finally, discusses and compares these results. 5.2 Method An experimental setup similar to that used in the ARF experiments of Chapter 4 was used in these homogeneous phantom experiments. The SonixMDP ul- trasound platform was used with the L14–5/38 ultrasound probe, which was clamped to a metal strut on an anti-vibration table (the setup was the same as shown in Fig. 4.1). Six phantoms of three different gelatine concentrations were constructed. The concentrations were 2, 3 and 4 wt%, and two were made of each concentration to test the repeatability of the imaging techniques. The phantoms were labelled A1, A2, B1, B2, C1 and C2, with types A, B and C rep- resenting 2, 3 and 4 wt%, respectively. All contained 4 wt% cellulose, giving 68 5.3. Axial Relaxation Imaging Results equivalent absorption properties and hence constant ARF across all three phan- toms. The values of K and K ′ (see Equations (2.1) and (2.11)) could therefore be assumed to be equal for each phantom, as is necessary for comparing rela- tive results from transfer function and parameter estimation in axial relaxation imaging. Before acquiring data, phantoms were left unrefrigerated and allowed to reach room temperature (18°C). In each experiment, a layer of distilled wa- ter was poured into the phantom cup, and the surface of the ultrasound probe was dipped into this layer, so that approximately 8 mm of water separated the surfaces of probe and phantom. Default sequence parameters were used in all acquisitions. In lateral relax- ation imaging, relaxation responses were acquired at six lateral locations spaced every 0.9144 mm from x1 = 1.829 mm to x6 = 6.401 mm. This corresponds to a shift of three probe elements between each location. Note that x = 0 corre- sponds to the pushing location at the centre of the probe. All RF data for all six phantoms were acquired within two hours, and were later processed offline. Displacement estimation was performed with a window size of 1 mm and a search radius of 100 µm. 5.3 Axial Relaxation Imaging Results This section gives results from the application of axial relaxation imaging to the homogeneous phantoms. See Sections 2.2.4, 2.3.3 and 2.3.4 for details of this imaging method. 5.3.1 Displacement Estimation The results of displacement estimation for the six phantoms are shown in Fig. 5.1. Each row gives the spatiotemporal displacement profiles for the two phan- toms of a given gelatine concentration. The focal-depth displacement relaxation responses for all phantoms are shown in Fig. 5.2. The relaxation responses were obtained by averaging the responses over a depth of 2 mm (two windows) cen- tred at the focal depth. 69 5.3. Axial Relaxation Imaging Results Figure 5.1: The axial relaxation imaging displacement estimation results for the six homogeneous phantoms. The spatiotemporal displacement profiles for the two phantoms of a given gelatine concentration phantoms are shown in each row. Time (ms) D ep th  (m m ) Displacement (µm) for phantom A1   0 5 10 15 20 25 10 15 20 25 30 35 0 5 10 15 20 25 30 35 40 45 50 Time (ms) D ep th  (m m ) Displacement (µm) for phantom A2   0 5 10 15 20 25 10 15 20 25 30 35 0 10 20 30 40 50 60  Time (ms) D ep th  (m m ) Displacement (µm) for phantom B1   0 5 10 15 20 25 10 15 20 25 30 35 0 5 10 15 20 Time (ms) D ep th  (m m ) Displacement (µm) for phantom B2   0 5 10 15 20 25 10 15 20 25 30 35 0 5 10 15 20  Time (ms) D ep th  (m m ) Displacement (µm) for phantom C1   0 5 10 15 20 25 10 15 20 25 30 35 -2 0 2 4 6 8 10 12 Time (ms) D ep th  (m m ) Displacement (µm) for phantom C2   0 5 10 15 20 25 10 15 20 25 30 35 -2 0 2 4 6 8 10 12  70 5.3. Axial Relaxation Imaging Results Figure 5.2: The focal-depth displacement relaxation responses obtained with axial relaxation imaging for the three homogeneous phantom types. 0 5 10 15 20 25 30 -10 0 10 20 30 40 50 60 Time (ms) D i s p l a c e m e n t  ( µ m ) Focal-depth displacement relaxation responses   Phantom A Phantom B Phantom C  5.3.2 Transfer Function Estimation The displacement estimation results of Section 5.3.1 were used with equation (2.4) to give the relative force-displacement transfer function estimates for each phantom. The Nyquist plots for these estimates are shown in Fig. 5.3 for all six phantoms and over the frequency range 0–1000 Hz. Equivalent magnitude and phase plots are shown in Fig. 5.4 from 0–800 Hz. Note that all results are normalised to the mean of the 0 Hz values of the type-A phantoms to give a clearer comparison. No information is lost, as, by (2.4), the transfer function estimates are already relative an unknown factor, K . 5.3.3 Parameter Estimation Applying the parameter estimation method of Section 2.3.4 to the axial relax- ation imaging transfer function estimation results of Section 5.3.2 gave esti- mates of relative values of the parameters of the described Voigt model, µr and ηr , the flow index of the power-law model, n, and the relaxation time constant 71 5.3. Axial Relaxation Imaging Results Figure 5.3: The Nyquist plots of the relative force-displacement transfer function estimates for the three phantoms types. Results are normalised to the mean of the 0 Hz values of the type-A phantoms. -0.2 0 0.2 0.4 0.6 0.8 1 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 Real I m a g i n a r y Nyquist plots of relative force-displacement transfer function   estimates  Phantom A Phantom B Phantom C 0 Hz0 Hz1000 Hz 300 Hz 100 Hz 100 Hz 1.2 0 Hz  Figure 5.4: The magnitude (left) and phase (right) plots of the relative force- displacement transfer function estimates for the three phantom types. Results are normalised to the mean of the 0 Hz values of the type-A phantoms. 102 10-1 100 Frequency (Hz) ative force-displacement transfer fun R e l a t i v e M a g n i t u d e Magnitude plots of rel ction estimates   0 Phantom A Phantom B Phantom C 102 -3 -2.5 -2 -1.5 -1 -0.5 Frequency (Hz) P h a s e  ( r a d ) Phase plots of relative force-displacement transfer function estimates   Phantom A Phantom B Phantom C  72 5.3. Axial Relaxation Imaging Results Figure 5.5: The relative elastic (left) and viscous (right) parameter estimation results for the three phantom types. Results are shown for both phantoms of each type. The elasticity results are normalised to the mean elasticity of the type- A phantoms. The viscosity results are normalised to the mean viscosity of the 100 Hz components of the type-A phantoms. Black bars indicate the relative Young’s moduli empirically predicted by (2.19). The power-law fits to the viscous parameter plots are shown as black dashed lines (right). A1 A2 B1 B2 C1 C2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Rel R e l a t i v e E l a s t i c i t y ative elasticity parameter estimates Phantom 0 100 200 300 400 500 600 700 800 0 0.5 1 1.5 2 2.5 3 3.5 4 Relative viscosity estimates with frequency Frequency (Hz) R e l a t i v e V i s c o s i t y   Phantom A Phantom B Phantom C Pow er Law  τ. The results for the relative elastic (µr ) and viscous (ηr ) parameters are shown in Fig. 5.5, together with the power-law fits to the viscous parameter plots. The elasticity results are normalised to the mean elasticity of the type-A phantoms, while the viscosity results are normalised to the mean viscosity of the 100 Hz components of the type-A phantoms. As for the transfer function estimation re- sults, this does not affect the outcome, due to the already relative nature of the estimates, but allows for better comparison. The relative Young’s moduli pre- dicted by the empirical formula, (2.19), are shown as black bars on the elasticity results. These values were also normalised for comparison. The power-law fits, plotted along with the viscous parameter plots in Fig. 5.5, were performed us- ing a linear least-squares fit to logarithmic plots of the relative viscous parame- ter with frequency. The fit was performed on the mean of the plots for the two phantoms of each type, and over the 100–400 Hz frequency range, where data appeared most reliable. The resulting flow index estimates were nA = −0.10, nB = −0.14 and nC = −0.18 for phantom types A, B and C, respectively. The re- laxation time constant results for the three phantom types are shown in Fig. 5.6. 73 5.4. Lateral Relaxation Imaging Results Figure 5.6: The relaxation time constant estimates with frequency for the three phantom types. 0 100 200 300 400 500 600 700 800 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 Time constant estimates with frequency Frequency (Hz) T i m e C o n s t a n t  ( m s )   Phantom A Phantom B Phantom C  5.4 Lateral Relaxation Imaging Results This section gives results from the application of lateral relaxation imaging to the homogeneous phantoms. See Sections 2.2.4, 2.3.5 and 2.3.6 for details of this imaging method. 5.4.1 Displacement Estimation After displacement estimation, the spatiotemporal displacement profiles for all phantoms and lateral imaging locations were used to provide focal-depth relax- ation responses by taking the mean of each profile over an axial range of 2 mm (two windows) centred on the focal depth. Each of these relaxation responses were filtered with the zero-phase FIR filter described in Section 2.3.5. The re- sulting plots for each phantom are shown in Fig. 5.7. 74 5.4. Lateral Relaxation Imaging Results Figure 5.7: The focal-depth displacement relaxation responses obtained with lateral relaxation imaging at six lateral locations for the three homogeneous phantom types, A (top), B (middle) and C (bottom). 30 0 5 10 15 20 25 30 -5 0 5 10 15 20 25 D i s p l a c e m e n t  ( µ m ) Lateral relaxation responses at focal depth for type-A phantoms   1.88 mm 2.82 mm 3.77 mm 4.71 mm 5.65 mm 6.59 mm Time (ms)  12 0 5 10 15 20 25 30 -2 0 2 4 6 8 10 Time (ms) D i s p l a c e m e n t  ( µ m ) Lateral relaxation responses at focal depth for type-B phantoms   1.88 mm 2.82 mm 3.77 mm 4.71 mm 5.65 mm 6.59 mm  8 0 5 10 15 20 25 30 -2 -1 0 1 2 3 4 5 6 7 Time (ms) D i s p l a c e m e n t  ( µ m ) Lateral relaxation responses at focal depth for type-C phantoms   1.88 mm 2.82 mm 3.77 mm 4.71 mm 5.65 mm 6.59 mm  75 5.4. Lateral Relaxation Imaging Results Figure 5.8: A plot of lateral location versus time of peak displacement for all six phantoms. Crosses show estimates obtained from lateral relaxation responses, while lines show linear least-squares fits for each phantom to the estimates. 1 2 3 4 5 6 7 x 10-3 0 1 2 3 4 5 6 7 8 9 x 10-3 Lateral Location (m) T i m e o f  p e a k d i s p l a c e m e n t  ( s ) Lateral location vs time of peak displacement   A1 A2 B1 B2 C1 C2  5.4.2 ShearWave Velocity and Parameter Estimation The lateral relaxation responses were used to obtain estimates of times of peak displacement at each lateral location for each phantom. The results are plotted in Fig. 5.8, along with the linear least-squares fit to the data for each phantom. The shear wave velocity estimates, cs , were obtained as the inverse of the gradients of the least-squares fits. The Young’s modulus, E , of each phantom was then calculated using (2.18). Relative Young’s modulus values, Er , were also calculated by dividing by the mean absolute value of the type-A phantoms. Ab- solute and relative Young’s modulus values predicted by the empirical formula (2.19), Eg and Eg ,r , respectively, were calculated for comparison. All these re- sults are given in Table 5.1, together with the relative elasticity parameter values from the axial relaxation imaging results of Section 5.3.3. 76 5.5. Discussion Table 5.1: The shear wave velocity (cs) and Young’s modulus (E) estimation re- sults of lateral relaxation imaging on the six homogeneous phantoms. The em- pirically predicted Young’s moduli are also shown (Eg ). Relative Young’s modu- lus or elastic parameter values from lateral relaxation imaging (Er ), the empir- ical formula (Eg ,r ) and axial relaxation imaging (µr ) are given in the last three columns for comparison. For each estimation method, relative values are nor- malised by the mean value of the type-A phantoms. Phantom cs (m/s) E (kPa) Eg (kPa) Er Eg ,r µr A1 0.70 1.46 1.78 1.01 1 1.03 A2 0.69 1.42 1.78 0.99 1 0.97 B1 1.04 3.24 4.16 2.25 2.34 2.35 B2 1.06 3.39 4.16 2.36 2.34 2.26 C1 1.48 6.61 7.58 4.60 4.26 4.63 C2 1.41 6.00 7.58 4.17 4.26 4.62 5.5 Discussion 5.5.1 Axial Relaxation Imaging The spatiotemporal displacement profiles of Fig. 5.1 show the expected peak displacement at the 20 mm focal depth and the subsequent relaxation be- haviour. It can be seen from the colour bars at the side of each profile that the peak displacement decreases as the gelatine concentration increases, showing that the stiffnesses of the phantoms increase with gelatine concentration. This can be seen more readily in the focal-depth relaxation responses of Fig. 5.2. Phantoms types A, B and C have peak displacements of approximately 56, 23 and 12 µm, respectively. Excellent consistency is observed between the two phantoms of each type in these relaxation responses and in all subsequent results. This indicates that the experimental method is not prone to any significant variation between ac- quisitions, and that any variations are due to differences between transmission media alone. The transfer function estimation results of Figs. 5.3 and 5.4 clearly show that different transfer functions are obtained for the three different phantom types. In the Nyquist plots, the value at 0 Hz and the overall radius of the curve de- 77 5.5. Discussion creases with increasing gelatine concentration. The shapes of all the Nyquist plots match the semi-circular curve of the Voigt model plot of Fig. 2.11, indi- cating that the Voigt model is a suitable choice for the parameter estimation approach used in this work. Note that the estimated Nyquist plots cross the imaginary axis at high frequencies, unlike the Voigt model. This indicates that though the Voigt model is suitable at lower frequencies, a higher-order model may be required to model the phantoms at higher frequencies. In all the magni- tude and phase plots, a low-pass response is observed, agreeing with the Voigt model. Characteristics that distinguish between phantoms as gelatine concen- trations decrease include a larger 0 Hz magnitude component, a steeper mag- nitude roll-off and lower cut-off frequency, and a steeper phase roll-off. This result also agrees with the Voigt model, as a medium with a lower elastic param- eter gives a magnitude response with a larger 0 Hz component, steeper roll-off and lower cut-off frequency, and a phase response with a steeper roll-off (see equations (2.9) and (2.10)). Although the viscous parameter of the Voigt model has the opposite effect, the elastic component has a greater dependence than the viscous on gelatine concentration, and is therefore the dominant parame- ter, as mentioned in Section 2.4.1. These results show that parameters derived directly from the transfer function estimates could be useful in viscoelastic char- acterisation of transmission media. The estimated relative elastic parameter of the Voigt model, µr , can be seen in the left of Fig. 5.5 to increase with gelatine concentration, with average val- ues of 1, 2.31 and 4.63 for phantom types A, B and C, respectively. The results are very consistent for both phantoms of each type. The black bars indicating the empirical Young’s modulus predictions of (2.19) – 1, 2.34 and 4.26 – are very well matched to the estimated parameters, giving a maximum error of 8.7%, for phantom C1. This shows the estimated relative elastic parameters follow the expected square-law relationship with gelatine concentration. It indicates the axial relaxation imaging method gives reliable estimates of the relative elastic Voigt parameter of homogeneous gelatine phantoms. The estimated relative viscous parameter, ηr , shown on the right of Fig. 5.5, also show a general increase with increasing gelatine concentration. As ex- pected, the increase is not as large as that of the relative elastic parameter, with, 78 5.5. Discussion for example, mean relative values at 100 Hz of 1, 1.80 and 2.97 for phantom types A, B and C, respectively. The results are consistent for both phantoms of each type, particularly over 100–400 Hz. The estimated viscosity parameters also de- crease with frequency, agreeing with the expectations of the power-law model. The flow index values, n, of the power-law model fit (–0.10, –0.14 and –0.18 for phantoms A, B and C, respectively) show that the viscosity parameters of these phantoms have a relatively weak frequency dependence, less than that of results previously obtained for 12 and 18 wt% gelatine phantoms, where n =−0.66 and −0.65, respectively [14]. The values also show that the dependence becomes stronger with gelatine concentration, a result that is observable in the viscosity plots of Fig. 5.5. The flow index values are all negative, as expected for gelatine media. The black dashed lines in Fig. 5.5 show that the fitted power-law mod- els match the estimation results well over 100–400 Hz for type-C phantoms and 100–500 Hz for type-A and -B phantoms. The relaxation time constant plots of Fig. 5.6 show that relaxation time con- stant generally decreases with gelatine concentration. They also show the same frequency-dependent trends as the relative viscous parameter plots. The mean values at 100 Hz for the three phantom types are 1.9, 1.4 and 1.2 ms. These val- ues match the observed relaxation responses well, as can be seen by plotting the observed relaxation responses together with exponential decay curves using these relaxation time constants (Fig. 5.9). Overall, the axial relaxation imaging results show that the method works well to give measures that can be used to differentiate the viscoelastic properties of different homogeneous transmission media. The estimated relative elastic and viscous parameters follow trends with gelatine concentration that agree with predictions and results from literature. This indicates that, together with the flow index and relaxation time constant parameters, they could be used suc- cessfully in characterisation of soft tissue. 5.5.2 Lateral Relaxation Imaging All the lateral relaxation responses of Fig. 5.7 clearly reveal the lateral propa- gation of the shear wave front. Responses obtained further from the pushing 79 5.5. Discussion Figure 5.9: A comparison of the observed relaxation responses and the relax- ation responses generated with the 100 Hz component of the relaxation time constant estimate. Generated responses are assumed to have an initial displace- ment equal to that of the observed responses. 0 5 10 15 20 25 30 -10 0 10 20 30 40 50 60 Time (ms) D i s p l a c e m e n t  ( µ m ) Experimental and time constant estimate relaxation res  ponses  Phantom A Phantom B Phantom C Time Constant Results  location have a peak displacement that occurs later in time and with a lower am- plitude, showing that the shear wave is attenuated by the medium as it travels. It can be seen that as gelatine concentration increases, the peaks are grouped close together (indicating a higher shear wave velocity) and experience higher attenuation with distance. This indicates the elasticity and viscosity increase with gelatine concentration, as for previous results. All responses are consistent over the two phantoms of each type, as with the axial relaxation imaging results. This consistency is also seen in the plots of lateral location versus time of peak displacement in Fig. 5.8. The lines showing the linear least squares fits for each phantom show excellent agreement to the data, with a maximum norm of residuals of 0.00022 for phantom A1. The gradients decrease with gelatine concentration, and give the shear wave velocity estimates in Table 5.1, which increase with gelatine concentration. Table 5.1 also shows the results of Young’s modulus parameter estimation for lateral relaxation imaging, with comparison to the empirical formula (2.19) and the relative elasticity parameter results of axial relaxation imaging. The ab- 80 5.5. Discussion solute values from (2.19) are somewhat higher than the estimated values. As mentioned in Section 2.4.1, absolute values can vary significantly depending on the construction and measurement methods. Also, the gelatine phantoms in the work from which (2.19) is taken [20] contained a small amount (0.4 wt%) of formaldehyde, a cross-linking agent that increases the final stiffness of the phantom. This may explain why the empirical predictions are higher than the values obtained from lateral relaxation imaging. The relative values, however, show good agreement for all phantoms between all three methods (empirical prediction, axial relaxation imaging and lateral relaxation imaging), with a max- imum error of 10.8% between any two methods. Finally, by assuming the elasticity parameter of the Voigt model is approxi- mately equal to the Young’s modulus for a medium (µ = E), absolute values of the viscosity parameter, η, can be obtained for each phantom using (2.15). Mean values at 100 Hz for phantom types A, B and C are 2.6, 4.7 and 7.7 Pa·s, respec- tively. The results of application of lateral relaxation imaging to homogeneous phantoms support the results of axial relaxation imaging and provide further insight into the behaviour of the these phantoms in response to ARF pushing pulses. Using this method with axial relaxation imaging allows absolute in- stead of relative viscoelastic parameters to be extracted, which could be of use in quantitative characterisation of soft tissue. 81 Chapter 6 Conclusion 6.1 Summary and Contributions This work presented a new ARF imaging technique called axial relaxation imag- ing that gives the ability to extract frequency-dependent viscoelastic informa- tion regarding soft tissue-like media and shows potential for use in diagnostic tissue classification. The theoretical basis of axial relaxation imaging, and the companion lateral relaxation imaging method, was detailed, along with their implementation on the SonixMDP ultrasound platform with Texo and a PC with MATLAB. In the system, a series of 100, medium-duration, focused pushing pulses, delivered at a high PRF, displaces the transmission medium. After the pushing pulses, short-duration imaging pulses track the resulting displacement relaxation response, either at the pushing location (in axial relaxation imaging) or at a lateral location some distance from the pushing location (in lateral relax- ation imaging). The digitised RF echo data are transferred to a PC for processing, through which viscoelastic properties of the medium are obtained. Patient and equipment safety issues related to the ARF imaging system were assessed. Patient safety was related to FDA guidelines for diagnostic ultrasound. Using acoustic water tank measurements, it was found that the acoustic output from the ARF imaging method met limits for mechanical index and spatial-peak, pulse-average intensity, but did not meet the limit for spatial-peak, time-average intensity during the pushing phase. It was noted, however, that the limit was only exceeded temporarily, and that, so long as the number of sequence rep- etitions per second was kept below ten, the method should be safe for in vivo patient use. A method for preventing damage to the ultrasound probe by limit- ing the number of sequence repetitions per acquisition, based on consideration of delivered energy, was also presented. 82 6.1. Summary and Contributions The ARF-generating and displacement tracking capabilities of the sys- tem were initially tested on a homogeneous gelatine phantom, the results of which gave valuable information regarding the abilities and limitations of the SonixMDP ultrasound platform in ARF imaging. It was found that a series of 100, 48-cycle, 6.7 MHz pushing pulses, delivered at a PRF of 20 kHz from an aperture of 32 elements, gave reasonable displacements in the phantom without over-loading the SonixMDP power supply. Also, the 100 pushing pulses were ob- served to create a near-saturated displacement, validating the assumption in ax- ial relaxation imaging transfer function estimation that the system is at a steady state before the pushing pulses cease. A small overall increase in ARF magnitude and resulting displacement may be obtained by increasing the pushing PRF to 60 kHz and using a centre frequency of 8 MHz. Axial and lateral relaxation imaging were then performed on a number of homogeneous gelatine phantoms with different gelatine concentrations. The results showed that axial relaxation imaging gave relative elastic and viscous parameters with very good separation between phantoms of different gelatine concentration; the 3 and 4 wt% phantoms were approximately 2.3 and 4.4 times more elastic and 1.8 and 3.0 times more viscous, respectively, than the 2 wt% phantoms. The elastic parameter results agreed within a 10.8% error to empiri- cal predictions and lateral relaxation imaging estimates of relative Young’s mod- ulus. It was also shown that the relaxation time constant parameter gave good separation between phantom types, without being dependent on the magni- tude of the ARF and while still giving frequency-dependent information. The frequency dependence of the viscous parameter for each of the phantoms was evaluated using the flow index parameter of the power-law relationship. This gave values between –0.18 and –0.10 for the three phantom concentrations, showing a relatively weak frequency dependence that increased with gelatine concentration. This work has made several important contributions, as detailed below. • Novel method of transfer function and frequency-dependent viscoelas- tic parameter estimation: A novel new method of estimating the force- displacement transfer function from a single displacement relaxation re- 83 6.1. Summary and Contributions sponse is introduced in this work. The transfer function is related in the frequency domain to a Voigt model with a frequency-dependent viscous parameter, to give elastic, viscous and time constant parameter estimates at a range of frequencies. The frequency dependence of the viscous pa- rameter is further modelled with a power law, giving a flow index param- eter estimate. Considering the range of information obtained, the acqui- sition method is very simple and fast. Comparative ARF imaging methods that provide frequency-dependent viscoelastic information either require long periods of ARF excitation [28, 72], give only frequency components at multiples of the pushing frequency [9], or employ more complex, fast ul- trasound imaging techniques [13]. The continuous frequency content of the relaxation response allows the viscoelastic parameters obtained with this method to be extracted from a 30 ms imaging phase following a sin- gle 5 ms pushing phase. Periodic, swept-frequency methods use continu- ous wave ARF, which may be applied for over 500 ms to obtain an equiv- alent frequency response (e.g. [28]). The shorter excitation and moni- toring times result in a faster acquisition time and limit the exposure of the medium to heating-related biological effects. The use in this work of transfer function estimation, and the extraction of time constant and flow index parameters, is unique in ARF imaging methods. • SonixMDP implementation: To the best of the author’s knowledge, this work represents the first implementation of an ARF imaging technique based around the SonixMDP ultrasound platform. This is important, as the open research architecture of the SonixMDP will allow other re- searchers to follow this work and develop their own ARF imaging tech- niques. Previous ARF work has been performed on platforms with closed architectures or on custom-built equipment, either requiring an exclusive agreement with an commercial ultrasound platform developer, or an over- head of work in the construction of a custom ultrasound system. The de- scription of the Texo sequence-generating software and parameters, and the assessment of the abilities and limitations of the SonixMDP, should be very beneficial to other researchers looking to develop ARF techniques 84 6.2. Future Work with this platform. • Frequency-dependent viscoelastic information regarding low- concentration gelatine phantoms: The testing of the techniques on gelatine phantoms gives useful information regarding elastic, viscous and frequency-dependent properties of low-concentration gelatine phantoms. The low gelatine concentration of these phantoms makes other mechanical testing techniques difficult, as they are easily damaged. There is therefore not a lot of information available in literature for such phantoms. Additionally, the frequency range over which parameter information is obtained by the axial relaxation imaging technique reveals information at frequencies typically not accessible to mechanical testing techniques. 6.2 FutureWork The ARF imaging techniques presented in this thesis are relatively new develop- ments, and as such, a lot of future work is required to bring them to the stage where they could be applied in a clinical setting to aid in tissue characterisation. Some such work is detailed below. • SonixMDP power supply: The SonixMDP currently experiences power supply voltage droop when transmitting long-duration pushing pulses. If stiffer media are to be imaged, this issue will need to be resolved, as it cur- rently limits the displacement that can be generated. External hardware may be required to provide additional power during pushing. • SonixMDP pulse limitations: As with the power supply voltage droop, the current pushing pulse parameters operate near the limit of what can be programmed with Texo. Longer pulse durations may be required to gener- ate detectable displacements in stiffer media. Again, this may be achieved with external hardware. Note that it would require re-evaluation of acous- tic output levels. 85 6.2. Future Work • Validation of viscous and power-law parameters: The results of viscous and power-law parameter estimation presented in this work have received little validation beyond qualitative assessment of relative values in rela- tion to other works in the literature. A reference method for measuring these properties would be beneficial. • Viscoelastic models: Models other than the Voigt model may prove more suitable for modelling the relaxation responses of soft tissue-like media. A fuller exploration of such models would bring more perspective on this topic to the work. • Further development of lateral relaxation imaging: The current lateral re- laxation imaging technique simply mimics other techniques from the lit- erature. Initial work has begun on developing transfer function estima- tion techniques with this method. It holds promise, because, as seen in the results of this work, absolute rather than relative parameters can be obtained, which would be beneficial for quantitative tissue characterisa- tion. • 2D axial relaxation imaging: More work is required to develop 2D axial re- laxation imaging methods for heterogeneous media. The spatial variabil- ity of the acoustic intensity and hence ARF profile, particularly with axial depth, needs to be considered in relation to relative parameter estimates. • Explore error sources in the system: Several potential sources of error exist in the current ARF imaging system. For example, displacement estimation does not currently account for induced lateral displacements. The finite sampling frequency of relaxation responses means the transfer function and parameter estimation methods also contain errors, and these need to be more fully explored. • In vivo application: The motivation for this work is to develop a system for tissue characterisation. Therefore the methods need eventually to be applied in vivo. Although the current implementation suffers from a lack of pushing ability, it may be applicable to softer tissues. Candidates for such experiments should be explored. 86 Bibliography [1] S. R. Aglyamov, S. Park, Y. A. Ilinskii, and S. Y. Emelianov, “Ultrasound imag- ing of soft tissue shear viscosity,” in IEEE Ultrasonics Symposium, Hon- olulu, HI, October 2003, pp. 937–940. [2] J. Bercoff, T. Mickaël, and M. Fink, “Supersonic shear imaging: A new tech- nique for soft tissue elasticity mapping,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 51, no. 4, pp. 396–409, April 2004. [3] J. Bercoff, T. Mickaël, M. Muller, and M. Fink, “The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Con- trol, vol. 51, no. 11, pp. 1523–1536, November 2004. [4] R. T. Beyer, Nonlinear Acoustics. Naval Ship Systems Command, Depart- ment of the Navy, 1974. [5] R. Bracewell, The Fourier Transformand Its Applications, 3rd ed. New York: McGraw-Hill, 2000, ch. 4, p. 63. [6] Canadian Cancer Society’s Steering Committee, “Canadian cancer statis- tics 2010,” Toronto, April 2010. [7] S. Catheline, J.-L. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelkaram, and J. Culioli, “Measurment of viscoelastic properties of homogeneous soft solid using transient elastography: An inverse problem approach,” Journal of the Acoustical Society of America, vol. 116, no. 6, pp. 3734–3741, Decem- ber 2004. [8] S. Chen, M. Fatemi, and J. F. Greenleaf, “Remote measurement of material properties from radiation force induced vibration of an embedded sphere,” 87 Bibliography Journal of the Acoustical Society of America, vol. 112, no. 3, pp. 884–889, September 2002. [9] S. Chen, M. W. Urban, C. Pislaru, R. Kinnick, Y. Zheng, A. Yao, and J. F. Greenleaf, “Shearwave dispersion ultrasound vibrometry (SDUV) for mea- suring tissue elasticity and viscosity,” IEEE Transactions onUltrasonics, Fer- roelectrics and Frequency Control, vol. 56, no. 1, pp. 55–62, January 2009. [10] D. A. Christensen, Ultrasonic Bioinstrumentation. New York, NY: Wiley, 1988. [11] J.-P. Corriou, Process Control: Theory and Applications. London: Springer, 2004, ch. 5, pp. 192–193. [12] P. G. M. de Jong, T. Arts, A. P. G. Hoeks, and R. S. Reneman, “Determination of tissue motion velocity by correlation interpolation of pulsed ultrasonic echo signals,” Ultrasonic Imaging, vol. 12, pp. 84–98, 1990. [13] T. Deffieux, G. Montaldo, M. Tanter, and M. Fink, “Shear wave spectroscopy for in vivo quantification of human soft tissues visco-elasticity,” IEEETrans- actions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 28, no. 3, pp. 313–322, March 2009. [14] H. Eskandari, S. E. Salcudean, and R. Rohling, “Viscoelastic parameter esti- mation based on spectral analysis,” IEEE Transactions on Ultrasonics, Fer- roelectrics, and Frequency Control, vol. 55, no. 7, pp. 1611–1625, July 2008. [15] M. Fatemi and J. F. Greenleaf, “Ultrasound-stimulated vibro-acoustic spec- trography,” Science, vol. 280, no. 5360, pp. 82–85, April 1998. [16] S. M. Frew, “ARF power supply testing,” University of British Columbia, Tech. Rep., May 2010. [17] L. Gao, K. J. Parker, R. M. Lerner, and S. F. Levinson, “Imaging of the elastic properties of tissue—a review,”Ultrasound inMedicine andBiology, vol. 22, no. 8, pp. 959–977, May 1996. 88 Bibliography [18] S. Girnyk, A. Barannik, E. Barannik, V. Tovstiak, A. Marusenko, and V. Volokhov, “The estimation of elasticity and viscosity of soft tissues in vitro using the data of remote acoustic palpation,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 211–219, February 2006. [19] B. B. Guzina, K. Tuleubekov, D. Liu, and E. S. Ebbini, “Viscoelastic charac- terization of thin tissues using acoustic radiation force and model-based inversion,” Ultrasound in Medicine and Biology, vol. 54, no. 13, pp. 4089– 4112, June 2009. [20] T. J. Hall, M. Bilgen, M. F. Insana, and T. A. Krouskop, “Phantom materials for elastography,” IEEE Transactions on Ultrasonics, Ferroelectrics and Fre- quency Control, vol. 44, no. 6, pp. 1355–1365, November 1997. [21] C. Joly-Duhamel, D. Hellio, and M. Djabourov, “All gelatin networks: 1. Bio- diversity and physical chemistry,” Langmuir, vol. 18, pp. 7208–7217, 2002. [22] M. Z. Kiss, T. Varghese, and T. J. Hall, “Viscoelastic characterization of in vitro canine tissue,” Physics inMedicine and Biology, vol. 49, pp. 4207–4218, 2004. [23] E. E. Konofagou and K. Hynynen, “Localized harmonic motion imaging: Theory, simulations and experiments,” Ultrasound in Medicine and Biol- ogy, vol. 29, no. 10, pp. 1405–1413, October 2003. [24] E. E. Konofagou, J. D’hooge, and J. Ophir, “Myocardial elastography—a fea- sibility study in vivo,” Ultrasound in Medicine and Biology, vol. 28, no. 4, pp. 475–482, April 2002. [25] T. A. Krouskop, T. M. Wheeler, F. Kallel, B. S. Garra, and T. J. Hall, “Elas- tic moduli of breast and prostate tissues under compression,” Ultrasonic Imaging, vol. 20, no. 4, pp. 260–274, October 1998. [26] S. A. Kruse, J. A. Smith, A. J. Lawrence, M. A. Dresner, A. Manduca, J. F. Greenleaf, and R. L. Ehman, “Tissue characterization using magnetic reso- nance elastography: Preliminary results,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1579–1590, June 2000. 89 Bibliography [27] C. P. Lee and T. G. Wang, “Acoustic radiation pressure,” Journal of the Acous- tical Society of America, vol. 94, no. 2, pp. 1099–1109, August 1993. [28] D. Liu and E. S. Ebbini, “Viscoelastic tissue property measurement using high frequency ultrasound,” in IEEE International Symposium on Biomed- ical Imaging, Arlington, VA, April 2007, pp. 892–895. [29] D. Liu and E. S. Ebbini, “Viscoelastic property measurement in thin tis- sue constructs using ultrasound,” IEEE Transactions on Ultrasonics, Ferro- electrics and Frequency Control, vol. 55, no. 2, pp. 368–383, February 2008. [30] E. L. Madsen, J. A. Zagzebski, R. A. Banjavie, and R. E. Jutila, “Tissue mim- icking materials for ultrasound phantoms,” Medical Physics, vol. 5, pp. 391– 394, 1978. [31] E. L. Madsen, J. A. Zagzebski, R. A. Banjavie, and M. M. Burlew, “Phantom material and method,” U.S. Patent 4,277,367, 1981. [32] E. L. Madsen, J. A. Zagzebski, and G. R. Frank, “Oil-in-gelatin dispersions for use as ultrasonically tissue-mimicking materials,” Ultrasound in Medicine and Biology, vol. 8, no. 3, pp. 277–287, October 1982. [33] E. L. Madsen, G. R. Frank, T. A. Krouskop, T. Varghese, F. Kallel, and J. Ophir, “Tissue-mimicking oil-in-gelatin dispersions for use in hetero- geneous elastography phantoms,” Ultrasonic Imaging, vol. 25, pp. 17–38, 2003. [34] R. G. Maev, R. E. Green, Jr., and A. M. Siddiolo, “Review of advanced acous- tical imaging techniques for nondestructive evaluation of art objects,” Re- search in Nondestructive Evaluation, vol. 17, pp. 191–204, 2006. [35] F. W. Mauldin Jr., O. B. Davis, M. A. Haider, E. G. Loboa, T. W. Pfeiler, and C. M. Gallippi, “On the potential of combined ARFI and elastography to improve differentiation of material structure in viscoelastic tissue,” in IEEE Ultrasonics Symposium, New York, NY, October 2007, pp. 2040–2045. 90 Bibliography [36] S. A. McAleavey, K. R. Nightingale, D. L. Stutz, S. J. Hsu, and G. E. Trahey, “Image reconstruction with acoustic radiation force induced shear waves,” in Proceedings of SPIE, vol. 5035, San Diego, CA, February 2003, pp. 223– 234. [37] M. Moradi, P. Mousavi, R. Rohling, and P. Abolmaesumi, “Tissue typing with ultrasound RF time series: Phantom studies,” in Proceedings of SPIE, vol. 7265, Lake Buena Vista, FL, February 2009, pp. 726 516-1–726 516-10. [38] L. A. Negron, F. Viola, E. P. Black, C. A. Toth, and W. F. Walker, “Development and characterization of a vitreous mimicking material for radiation force imaging,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 49, no. 11, pp. 1543–1551, November 2002. [39] K. R. Nightingale, P. J. Kornguth, W. F. Walker, B. A. McDermott, and G. E. Trahey, “A novel ultrasonic technique for differentiating cysts from solid le- sions: Preliminary results in the breast,” Ultrasound in Medicine and Biol- ogy, vol. 21, no. 6, pp. 745–751, 1995. [40] K. R. Nightingale, M. L. Palmeri, R. W. Nightingale, and G. E. Trahey, “On the feasibility of remote palpation using acoustic radiation force,” Journal of the Acoustical Society of America, vol. 110, no. 1, pp. 625–634, July 2001. [41] K. R. Nightingale, R. Bentley, and G. E. Trahey, “Observations of tissue re- sponse to acoustic radiation force: Opportunities for imaging,” Ultrasonic Imaging, vol. 24, no. 3, pp. 129–138, July 2002. [42] K. R. Nightingale, M. S. Soo, R. W. Nightingale, and G. E. Trahey, “Acoustic radiation force impulse imaging: In vivo demonstration of clinical feasibil- ity,” Ultrasound in Medicine and Biology, vol. 28, no. 2, pp. 227–235, Febru- ary 2002. [43] K. R. Nightingale, S. A. McAleavey, and G. E. Trahey, “Shear-wave generation using acoustic radiation force: In vivo and ex vivo results,” Ultrasound in Medicine and Biology, vol. 29, no. 12, pp. 1715–1723, December 2003. 91 Bibliography [44] K. R. Nightingale, L. Zhai, J. J. Dahl, K. D. Frinkley, and M. L. Palmeri, “Shear wave velocity estimation using acoustic radiation force impulsive excita- tion in liver in vivo,” in IEEE Ultrasonics Symposium, Vancouver, BC, Octo- ber 2006, pp. 1156–1160. [45] W. L. Nyborg, “Acoustic streaming,” in Physical Acoustics, W. P. Mason, Ed. New York: Academic Press, 1965, vol. IIB, pp. 265–331. [46] J. Ophir, I. Céspedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: A quantitative method for imaging the elasticity of biological tissues,” Ultra- sonic Imaging, vol. 13, no. 2, pp. 111–134, April 1991. [47] J. Ophir, B. S. Garra, F. Kallel, E. E. Konofagou, T. A. Krouskop, R. Righetti, and T. Varghese, “Elastographic imaging,” Ultrasound in Medicine and Bi- ology, vol. 26, Supplement 1, pp. S23–S29, 2000. [48] M. Orescanin, K. S. Toohey, and M. F. Insana, “Material properties from acoustic radiation force step response,” Journal of the Acoustical Society of America, vol. 125, no. 5, pp. 2928–2936, May 2009. [49] M. L. Palmeri and K. R. Nightingale, “On the thermal effects associated with radiation force imaging of soft tissue,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 51, no. 5, pp. 551–565, May 2004. [50] M. L. Palmeri, A. C. Sharma, R. R. Bouchard, R. W. Nightingale, and K. R. Nightingale, “A finite-element method model of soft tissue response to im- pulsive acoustic radiation force,” IEEE Transactions on Ultrasonics, Ferro- electrics and Frequency Control, vol. 52, no. 10, pp. 1699–1712, October 2005. [51] M. L. Palmeri, S. A. McAleavey, K. L. Fong, G. E. Trahey, and K. R. Nightin- gale, “Dynamic mechanical response of elastic spherical inclusions to im- pulsive acoustic radiation force excitation,” IEEE Transactions on Ultra- sonics, Ferroelectrics and Frequency Control, vol. 53, no. 11, pp. 2065–2079, November 2006. 92 Bibliography [52] M. L. Palmeri, M. H. Wang, J. J. Dahl, K. D. Frinkley, and K. R. Nightingale, “Quantifying hepatic shear modulus in vivo using acoustic radiation force,” Ultrasound in Medicine and Biology, vol. 34, no. 4, pp. 546–558, April 2008. [53] K. J. Parker, L. S. Taylor, S. Gracewski, and D. J. Rubens, “A unified view of imaging the elastic properties of tissue,” Journal of the Acoustical Society of America, vol. 117, no. 5, pp. 2705–2712, January 2005. [54] C. J. Pavlin, M. D. Sherar, and F. S. Foster, “Subsurface ultrasound micro- scopic imaging of the intact eye,” Ophthalmology, vol. 97, no. 2, pp. 244– 250, February 1990. [55] L. Rayleigh, “On the pressure of vibrations,” Philosophical Magazine, vol. 3, no. 15, pp. 338–346, March 1902. [56] L. Rayleigh, “On the momentum and pressure of gaseous vibrations, and on the connexion with the virial theorem,” Philosophical Magazine, vol. 10, no. 57, pp. 364–374, September 1905. [57] D. W. Rickey, P. A. Picot, D. A. Christopher, and A. Fenster, “A wall-less ves- sel phantom for doppler ultrasound studies,” Ultrasound in Medicine and Biology, vol. 21, no. 9, pp. 1163–1176, 1995. [58] O. V. Rudenko, A. P. Sarvazyan, and S. Y. Emelianov, “Acoustic radiation force and streaming induced by focused nonlinear ultrasound in a dissipative medium,” Journal of the Acoustical Society of America, vol. 99, no. 5, pp. 2791–2798, May 1996. [59] A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov, “Shear wave elasticity imaging: A new ultrasonic technology of medical diagnostics,” Ultrasound in Medicine and Biology, vol. 24, no. 9, pp. 1419–1435, November 1998. [60] R. Seip and E. S. Ebbini, “Noninvasive estimation of tissue temperature re- sponse to heating fields using diagnostic ultrasound,” IEEE Transactions on Biomedical Engineering, vol. 42, no. 8, pp. 828–839, August 1995. 93 Bibliography [61] R. Seip, P. VanBaren, C. A. Cain, and E. S. Ebbini, “Noninvasive real-time multipoint temperature control for ultrasound phased array treatments,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 43, no. 6, pp. 1063–1073, November 1996. [62] R. Sinkus, K. Siegmann, T. Xydeas, M. Tanter, C. Claussen, and M. Fink, “MR elastography of breast lesions: Understanding the solid/liquid dual- ity can improve the specificity of contrast-enhanced MR mammography,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1135–1144, December 2007. [63] T. Sugimoto, S. Ueha, and K. Itoh, “Tissue hardness measurement using the radiation force of focused ultrasound,” in IEEE Ultrasonics Symposium, Honolulu, HI, December 1990, pp. 1377–1380. [64] T. L. Szabo, Diagnostic Ultrasound Imaging: Inside Out. Burlington, MA: Elsevier Academic Press, 2004. [65] M. Tanter, J. Bercoff, A. Athanasiou, T. Deffieux, J.-L. Gennisson, G. Mon- taldo, M. Muller, A. Tardivon, and M. Fink, “Quantitative assessment of breast lesion viscoelasticity: Initial clinical results using supersonic shear imaging,”Ultrasound inMedicine andBiology, vol. 34, no. 9, pp. 1373–1386, September 2008. [66] Ultrasonix Medical Corporation, “Transducer specification sheet,” Ultra- sonix Medical Corp., August 2009. [67] Ultrasonix Medical Corporation, “Texo documentation 2.0,” Ultrasonix Medical Corp., July 2010. [Online]. Available: wikisonix/api/texo/2.0/ [68] M. W. Urban, S. Chen, and J. F. Greenleaf, “Error in estimates of tissue ma- terial properties from shear wave dispersion ultrasound vibrometry,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 56, no. 4, pp. 748–758, April 2009. 94 Bibliography [69] M. W. Urban, G. T. Silva, M. Fatemi, and J. F. Greenleaf, “Multifrequency vibro-acoustography,” IEEE Transactions on Medical Imaging, vol. 25, no. 10, pp. 1284–1295, October 2006. [70] M. W. Urban, M. Fatemi, and J. F. Greenleaf, “Modulation of ultrasound to produce multifrequency radiation force,” Journal of the Acoustical Society of America, vol. 127, no. 3, pp. 1228–1238, March 2009. [71] U.S. Food and Drug Administration, Information for Manufacturers Seeking Marketing Clearance of Diagnostic Ultrasound Systems and Trans- ducers, U.S. Food and Drug Administration Std., September 2008. [Online]. Available: DeviceRegulationandGuidance/GuidanceDocuments/UCM070911.pdf [72] J. Vappou, C. Maleke, and E. E. Konofagou, “Quantitative viscoelastic pa- rameters measured by harmonic motion imaging,” Physics inMedicine and Biology, vol. 54, no. 11, pp. 3579–3594, June 2009. [73] F. Viola and W. F. Walker, “Imaging viscoelastic properties of the vitreous,” in IEEE Ultrasonics Symposium, Atlanta, GA, October 2001, pp. 1623–1626. [74] F. Viola, F. W. Mauldin Jr., S. P. Tropello, B. G. Macik, M. B. Lawrence, and W. F. Walker, “Sonorheometry: A new method for assessing coagulation po- tential,” in IEEE Ultrasonics Symposium, New York, NY, October 2007, pp. 1001–1004. [75] W. F. Walker, R. E. Davidsen, and C. A. Toth, “Application of acoustic ra- diation force in ophthalmic ultrasound,” in IEEE Ultrasonics Symposium, Toronto, ON, October 1997, pp. 1291–1295. [76] M. H. Wang, M. L. Palmeri, N. C. Rouze, K. R. Nightingale, and M. A. Hob- son, “Investigating the effects of viscosity on focused, impulsive, acoustic radiation force induced shear wave morphology,” in IEEE Ultrasonics Sym- posium, Beijing, China, November 2008, pp. 647–650. [77] P. Webb, I. Gibson, and C. Wykes, “Robot guidance using ultrasonic arrays,” Journal of Robotic Systems, vol. 11, no. 8, pp. 681–692, 1994. 95 [78] T. A. Whittingham, “Tissue harmonic imaging,” European Radiology, vol. 9, Supplement 3, pp. S323–S326, 1999. [79] Y. Yamakoshi, J. Sato, and T. Sato, “Ultrasonic imaging of internal vibration of soft tissue under forced vibration,” IEEE Transactions onUltrasonics, Fer- roelectrics and Frequency Control, vol. 17, no. 2, pp. 45–53, March 1990. [80] R. Zahiri-Azar, “Acoustic radiation force on SonixRP,” Ultrasonix Medical Corp., Tech. Rep., April 2009. [81] L. Zhai, “Imaging and characterizing human prostates using acoustic radi- ation force,” Ph.D, Duke University, Durham, NC, 2009. [82] L. Zhai, S. J. Hsu, R. R. Bouchard, and K. R. Nightingale, “A combined ARFI sequence for 2D displacement imaging and shear wave velocity mapping,” in IEEE Ultrasonics Symposium, Beijing, China, November 2008, pp. 2013– 2016. [83] L. Zhai, J. Madden, W. C. Foo, M. L. Palmeri, V. Mouraviev, T. J. Polascik, and K. R. Nightingale, “Acoustic radiation force impulse imaging of human prostates ex vivo,” Ultrasound in Medicine and Biology, vol. 36, no. 4, pp. 576–588, April 2010. 96 Appendix A Derivation of an Expression for the Acoustic Radiation Force This appendix provides a derivation of the equation for the ARF per unit vol- ume, as given in Equation (1.1) of Chapter 1. The derivation follows those given by Nyborg [45] and Beyer [4] for acoustic streaming. These works derive gen- eral expressions for the ARF, then simplify for the case of a plane wave. Here, however, the plane wave simplification is made at the beginning. Firstly, assume a plane ultrasonic wave travelling in a viscous liquid medium which is unbounded, homogeneous and isotropic. The motion is governed by the Navier-Stokes equation for conservation of momentum, which, for the as- sumed form of the wave and medium, is given by ρ ( ∂v ∂t + v ∂v ∂z ) =−∂p ∂z + (η′+ 4 3 η) ∂2v ∂z2 . (A.1) Here, ρ is the density of the medium, v is velocity, t is time, z is the spatial co- ordinate for the direction of propagation (assumed to be axial), p is pressure, and η′ and η are the shear and bulk viscosity coefficients, respectively (note that these are not the same as the viscous parameter, η, used throughout the rest of this thesis). The terms on the left-hand side of (A.1) are related to inertia while those on the right-hand side represent the net force per unit volume due to stresses: FS =−∂p ∂z + (η′+ 4 3 η) ∂2v ∂z2 . (A.2) The equation of continuity, ∂ρ ∂t + ∂ ( ρv ) ∂z = 0, (A.3) 97 Appendix A. Derivation of an Expression for the Acoustic Radiation Force can be multiplied through by v and added to the left-hand side of (A.1) to give, after using the chain rule and simplifying, ∂ ( ρv ) ∂t +ρv ∂v ∂z + v ∂ ( ρv ) ∂z =−∂p ∂z + (η′+ 4 3 η) ∂2v ∂z2 . (A.4) Here, the last two terms of the left-hand side represent additional inertial com- ponents resulting from the fact that the medium has not been assumed to be incompressible. Together, they constitute an excess force per unit volume term, −FE , such that −FE = ρv ∂v ∂z + v ∂ ( ρv ) ∂z . (A.5) Note that this is in agreement with Newton’s second law, with the rate of change of the linear momentum per unit volume being equal to the net force per unit volume. I.e. from (A.2), (A.4) and (A.5), ∂ ( ρv ) ∂t = FS +FE . (A.6) Due to the nonlinear and complex nature of (A.4), perturbation theory is used to approach a solution. Here, the quantities p, ρ and v are expanded as a Taylor series about zero in the perturbation parameter, ², giving a sum of suc- cessively higher-order approximations. Thus, p = p0+²p1+²2p2+ . . . ρ = ρ0+²ρ1+²2ρ2+ . . . (A.7) v = ²v1+²2v2+ . . . , where the subscripts denote the order of the approximations. The zeroth-order values, p0 and ρ0, are the static values of pressure and density (while v0 = 0). Using zeroth-order (²0) terms in (A.4) gives ∂p0∂z = 0, such that p0 is a constant. Similarly, using the continuity equation (A.3), gives ∂ρ0∂t = 0, so that ρ0 is a con- stant. Considering only first-order (²1) terms in (A.4) gives 98 Appendix A. Derivation of an Expression for the Acoustic Radiation Force ρ0 ∂v1 ∂t =−∂p1 ∂z + (η′+ 4 3 η) ∂2v1 ∂z2 . (A.8) This is equivalent to the equation for an incompressible medium (the−FE terms disappear), and describes the linear wave propagation at the fundamental har- monic frequency. The velocity here is assumed to have the form of a damped harmonic travelling wave: v1 = vpe−αz sin(ωt −kz) , (A.9) where vp is the amplitude,ω is the fundamental angular frequency, k is the wave number, and α is the absorption coefficient. The absorption coefficient is given by α= ( η′+ 4 3 η ) ω2 2ρ0c3 , (A.10) where c = √ p0 ρ0 is the speed of sound. Now considering only the second-order (²2) terms, (A.4) becomes, after sim- plification, ∂ ∂t ( ρ0v2+ρ1v1 )+2ρ0v1 ∂v1 ∂z = −∂p2 ∂z + (η′+ 4 3 η) ∂2v2 ∂z2 . (A.11) The second-order quantities consist of time-independent and 2ω components, but only the time-independent components are of interest here. These can be obtained by averaging (A.11) in time over an integral number of cycles. This operation is represented by the 〈·〉 operator. The first term on the left-hand side becomes zero (see [4]): 〈 ∂ ∂t ( ρ0v2+ρ1v1 )〉= 0. (A.12) Hence, (A.11) becomes〈 2ρ0v1 ∂v1 ∂z 〉 = 〈−∂p2 ∂z + (η′+ 4 3 η) ∂2v2 ∂z2 〉 , (A.13) or, 〈−FE2〉 = 〈−∂p2 ∂z + (η′+ 4 3 η) ∂2v2 ∂z2 〉 , (A.14) 99 Appendix A. Derivation of an Expression for the Acoustic Radiation Force where 〈−FE2〉 is the time-average of the second-order terms of the excess force per unit volume, FE . The time-average excess force per unit volume is thus 〈FE2〉 = F =−2ρ0 〈 v1 ∂v1 ∂z 〉 , (A.15) where F is the desired body force due to ARF. It can be seen that this depends only on the first-order velocity. Substituting (A.9) gives F = 2ρ0v2pe−2αz 〈 αsin2 (ωt −kz)+k sin(ωt −kz)cos(ωt −kz)〉 = αρ0v2pe−2αz , (A.16) because the time-average of the sine-squared term isα/2, while that of the sine- cosine term is zero. This is the result obtained in Nyborg [45] and Beyer [4]. This expression can further be related to the time-average acoustic intensity (power per unit area) of the ultrasonic wave, 〈I 〉, using the relations I = pv and p = ρ0cv [10]. Thus, 〈I 〉 = 〈ρ0cv21〉 = ρ0cv2pe −2αz 2 , (A.17) and, substituting this into (A.16), F = 2α〈I 〉 c . (A.18) This is equation (1.1), as desired. 100 Appendix B The Relationship between the Step and Relaxation Responses The assumption that the relaxation response of a linear, time-invariant system can be converted to a step response, as used for equation (2.2), is justified here. Essentially, this consists of proving that K s (t ) = r (0)− r (t ), where K is a con- stant, s (t ) is the unit step response and r (t ) is the relaxation response. If s (t ) is the unit step response of a linear, time-invariant system, h (t ), to the Heaviside step function, θ (t ), then K s (t )=Kθ (t )⊗h (t ) , (B.1) where K is a constant and ⊗ denotes time-domain convolution. The relaxation response, r (t ), to a scaled, zero-going step function, K [1−θ (t )], input to h (t ), is r (t ) = K [1⊗h (t )−θ (t )⊗h (t )] = K ∞̂ −∞ h (t )dt −Kθ (t )⊗h (t ) = KH (ω= 0)−K s (t ) , (B.2) where H (ω= 0) is the transfer function of the system (the Fourier transform of h (t )) evaluated at the angular frequency ω= 0, and is hence a constant value. Now, from (B.2), r (0)=KH (ω= 0)−K s (0) at t = 0. The Heaviside step func- 101 Appendix B. The Relationship between the Step and Relaxation Responses tion has the property [5], θ (t )⊗ f (t )= tˆ −∞ f ( t ′ ) dt ′, (B.3) for a function f (t ). Using (B.1) and (B.3), gives K s (0) = K 0ˆ −∞ h ( t ′ ) dt ′ = 0, (B.4) assuming h (t ) = 0 for t ≤ 0. Thus, assuming the system is causal and does not have an instantaneous response, r (0) = KH (ω= 0). Then, from (B.2), r (t ) = r (0)−K s (t ) and K s (t )= r (0)− r (t ), as required. 102


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items