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.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
24-ubc_2011_spring_frew_samuel.pdf [ 5.71MB ]
Metadata
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
24-1.0071534-fulltext.txt
Citation
24-1.0071534.ris

Full Text

Viscoelastic Imaging Methods 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 imaging, 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 relative 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 characterisation.  ii  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table of Contents  ii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii  List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii  Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dedication  x  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 4 ARF Experiments  . . . . . . . . . . . . . . . . . . . . . . 55  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 successes 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, especially 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 images 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 relatively 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 scattering 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 characterisation 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 effectiveness 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 region 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 converting the movement into information that can aid in tissue characterisation. Techniques are typically grouped by the method of generating movement. Originally, 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 propagation 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 received 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 imaging [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 converted 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 provide for improving tissue characterisation. For example, cancerous breast tissue 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 provide 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 ultrasound transmission, the resulting movement occurs in a relatively localised region (roughly equivalent to the dimensions of the focal region of the transmission). This localised nature of ARF imaging means that boundary conditions can often be ignored, and that select regions can be palpated independently. 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 displacement 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 ultrasound 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 extract elastic, viscous and frequency-dependent parameters from soft ultrasonic media. The remainder of this chapter contains general background on ultrasonic 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 software, and the construction of the gelatine tissue phantoms used as the ultrasonic transmission media in the subsequent experiments. In Chapter 3, issues related to the safety of the methods to the potential patient and to the ultrasound equipment are studied. Chapter 4 describes initial ARF experiments performed 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 tissue 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 converts voltage signals produced by the ultrasound platform into ultrasonic pressure signals, typically using piezoelectric transducers. In pulse-echo ultrasound, short-duration (order of microseconds), sinusoidal pressure pulses are generated. The surface of the probe is placed in contact with the ultrasonic transmission medium, so that, with the aid of acoustic matching layers, the pressure pulse conducts into the medium. As the pressure wave of the pulse propagates 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 voltage signals, which are transferred to the ultrasound platform for processing and display. Essentially, the instantaneous amplitude of the received signal determines 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 information 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 transmitted pulse is delivered. A common way to shift this lateral location is by using a phased-array ultrasound probe, which contains numerous individual piezoelectric transducer elements, any or several of which can be electronically selected by the ultrasound platform. By using multiple elements in each pulse transmission 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. Lateral  Elevational  Axial  Display  Ultrasound Probe  Ultrasound Platform Pulse  Transmission Medium  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 techniques 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 absorbing object, the excess pressure results in a net force, the ARF. Further reading 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 media 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 derivation 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 propagates. Because higher density results in higher wave speed, positive compressional 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 direction 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 literature 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  Acoustic Radiation Force Impulse (ARFI)  Shear Wave Elasticity Imaging (SWEI)*  Nightingale et al . 2001  Sarvazyan et al., 1998  * indicates method uses monitoring of shear waves  Periodic  Vibro-acoustography Fatemi & Greenleaf , 1998  Acoustic Radiation Force Impulse (ARFI)*  Shear Wave Dispersion Vibrometry (SDUV)*  McAleavy et al., 2003  Chen et al., 2009  Supersonic Shear Imaging (SSI)*  Harmonic Motion Imaging (HMI)  Bercoff et al., 2004  Konofagou & Hynynen, 2003  Sonorheometry  High Frequency Harmonic Imaging  Walker et al., 1997  Liu & Ebbini 2007  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 resulting displacement or strain is measured after the system has reached a stable 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) imaging, developed by researchers at Duke University [40, 41, 42]. Here, several longduration (~20–60 µs), focused pushing pulses are delivered into the tissue, creating an “impulse” in ARF. In the static version of ARFI, the displacement response at the pushing location is monitored using standard pulse-echo ultrasound (from the same probe) during and after the pushing pulses, giving an approximate 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 regions 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 regions of the prostate anatomy, and detect changes due to prostate cancer, benign 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 frequencies, 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; however, 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 easily 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 positions. 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 propagation, with research by the same Duke University group undertaken to develop 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 methods 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 several 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 reconstruct 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 shortduration (< 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 sufficient 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 theoretical 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 parameters to be obtained. If multiple frequencies are used, frequency-dependencies can also be obtained. Like transient ARF techniques, methods may monitor displacements at the pushing location or at a lateral location away from the pushing 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 centre frequencies) interfere to produce an ARF that varies sinusoidally at the difference 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 complete frequency response over the swept range. This can then be related to theoretical 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 parameters. A similar periodic method is harmonic motion imaging [23], and recent 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 ultrasound, 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 frequencydependent 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 relaxation imaging and lateral relaxation imaging, use transient excitation methods. In axial relaxation imaging, in contrast to previous excitation methods, pushing 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 imaging pulses, the relaxation response after removal of ARF is monitored. The relaxation response, being a transient signal, has a frequency spectrum containing information about the medium over a range of frequencies. It is analysed in the frequency domain, applying transfer function techniques, and the frequency 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 single relaxation response (acquired in less than 30 ms). Other methods that give this frequency-dependent information typically rely on swept-frequency periodic ARF with longer excitation and monitoring times. Lateral relaxation imaging, 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. Particulars 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, occur in these blocks. Subsequent processing of data into meaningful information regarding the viscoelastic properties of the transmission medium and the display 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 relates 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 captured, 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 gen15  2.1. Ultrasound Platform and Probe  Figure 2.1: The overall system used in ARF imaging in this work. Texo Sequence Generating Software  Computing Device Data Import  Displacement Estimation  Ultrasound Platform  Axial Relaxation Imaging:  Lateral Relaxation Imaging:  Transfer Function Estimation  Shear Wave Velocity Estimation  Parameter Estimation  Display Conversion  Ultrasound Probe  Display  Transmission Medium  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 electrical pulses to the connected ultrasound probe and receives and processes the returned RF echoes. The processing includes filtering, sampling and quantisation 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 hardware 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 functions and parameters control the generation of ultrasound pulse sequences and processing of the received RF echo signals performed by the SonixMDP. Parameters 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 provided 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, either 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 parameters and image formation parameters. These are discussed further in Section 2.2.3. After setting the desired sequence parameters, the user can run the sequence, which begins the transmission of pulses and the sampling of the resulting 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. Mi cycles (tPDi)  Mp cycles (tPDp)  Imaging pulse  VTX  Pushing pulse  T c =1/f c  i.1  i.2  . . .  i.Ni  T PRi =1/f PRi Imaging phase  p.1  p.2  . . .  p.Np  i.1  . . .  T PRp =1/f PRp Pushing phase  Sequence loop  the “Stop sequence” command) or the callback function designed to prevent excessive 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 imaging phase, with the imaging pulse parameters, while the second part shows the pushing phase, with the pushing pulse parameters. The top part of the figure shows individual pulse characteristics, while the lower part shows sequence characteristics. The number of cycles in a pulse, M i or M p , 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 Pushing Pulse No. cycles PRF Transmit voltage Centre frequency No. pulses per sequence loop Imaging Pulse No. cycles PRF Transmit voltage Centre frequency No. pulses per sequence loop Image Formation Focal depth Acquisition depth Pushing transmit aperture Imaging transmit aperture Imaging receive aperture Pushing centre element Imaging centre element TGC top gain TGC bottom gain TGC bottom depth TGC middle gain TGC middle depth Digital gain Receive sampling frequency  Symbol  Default Value  Modifiable Range  Mp f P Rp VT X fc Np  48 20 kHz 47 V 6.7 MHz 100  1–48 1–200 kHz N/A 4–20 MHz 1–1000  Mi f P Ri VT X fc Ni  1 10 kHz 47 V 6.7 MHz 300  1–48 1–20 kHz N/A 4–20 MHz 1–1000  Df Da AT X p AT X i AR X i Xp Xi Gt Gb zb Gm zm GD fs  20 mm 40 mm 32 elements 32 elements 32 elements 64 64 40% 60% Df 50% 0.5D f 2500 20 MHz  10–50 mm 10–50 mm 0–128 0–128 0–64 0–128 0–128 0–100% 0–100% 10–50 mm 0–100% 10–50 mm 0–4095 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 ; acquisition depth, D a ; transmit/receive aperture, A T X /R X ; and imaging/pushing centre element, X i /p .  Lateral Location 0  1  Axial Depth  0  z  ATX/RX  x Ultrasound Probe  Xi/p 128  ...  ...  Probe Elements  Df Focus  Da  Transmission Medium  21  2.2. Texo Sequence Generating Software nusoidal cycles at the centre frequency, f c , making up the pulse duration, t P Di or t P D p . These parameters are related by t P D = M / f c , where the appropriate subscript should be applied to t P D 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, f P Ri or f P Rp , which is the inverse of the pulse repetition period (PRP), TP Ri or TP Rp , and by the number of pulses per sequence loop, Ni or N p (labelled in the figure as i.1 to i.Ni for imaging pulses and p.1 to p.N p 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, D a 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 transmission 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 media 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 A T X i , A T X p and A R X i (shown collectively as A T X /R X 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, X i and X p 22  2.2. Texo Sequence Generating Software (shown as X i /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 signals. The digital gain, G D , 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 frequency, f s , 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) , G t ; the gain and position at the maximum depth (the “bottom” of the curve), G b and z b , respectively; and the gain and depth of a mid-point between the top and bottom, G m and z m , 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 z b is set to the value of the focal depth, D f , while the default value of z m is half the value of z b . Note that depths beyond the maximum depth, z b , have the same gain, G b , 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 Section 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, G t , bottom gain, G b , bottom depth, z b , middle gain, G m , and middle depth, z m .  Gain (%)  100  zm  50  zb  Gm  Gt 0  10  Gb  20  30  40  Depth (mm)  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 lateral relaxation imaging comes in the imaging phase. In axial relaxation imaging, the lower-intensity imaging pulses monitor the axial movement at the same lateral 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 viscoelastic 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 ultrasound 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. Pushing  Axial Relaxation Imaging  Lateral Relaxation Imaging  x z  Ultrasound Probe  Transmission Medium  Push creates movement  Pushing location  Movement imaged at same location  Movement imaged at different 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 after the next pushing phase. This method, however, can cause issues if the movements 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, X i , is adjusted for each acquisition. Each output spatiotemporal 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 profiles. 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 re25  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 package 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 processing 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 estimation (for lateral relaxation imaging), parameter estimation and display conversion. 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 1 The 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. 4  1  RF echo amplitude with depth  x 10  0.8 0.6  Amplitude  0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1  0  5  10  15 20 25 Axial Depth (mm)  30  35  40  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 snapshot 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 estimate 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 func2 The 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, f P Ri . The individual RF echo signals can be up-sampled from the receive sample frequency, f s , before 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 normalised 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 u A , u B 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, giving 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 transmission medium at different times. B  C  A  Current RF echo  uA  uB  uC  B  C  A  18.5  Reference RF echo  19  19.5  20  20.5  z  Axial Depth (mm)  29  2.3. Data Processing Software  2.3.3 Axial Relaxation Imaging: Transfer Function Estimation This section describes the transfer function estimation that follows displacement 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, equivalently, 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 displacement, 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 scaling 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 discussed 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 pushing phase, which seems a reasonable assumption given that the pushing pulses stop and the imaging pulses that begin are of very low intensity in comparison. 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, s r (t ), is calculated as s r (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.  Force f(t)  Transfer Function H(ω)  Displacement u(t)  cause the derivative of a step function is an impulse function, the assumed linearity of the system allows the relative impulse response, h r (t ), to be obtained from the relative step response by differentiation in time. Thus, h r (t ) = K h (t ) =  d d s r (t ) = − r (t ) , dt dt  (2.3)  because r (0) is a constant. Here, h (t ) is the (non-scaled) impulse response. Finally, taking the Fourier transform of (2.3) gives the relative, force-displacement transfer function estimate, Hr (ω) = K H (ω) = K  U (ω) = F {h r (t )} , F (ω)  (2.4)  where ω is angular frequency, U (ω) and F (ω) are the Fourier transforms of u (t ) and f (t ), respectively, and F {·} represents the Fourier transform operation. The diagram of Fig. 2.9 shows the system diagram used for identifying the forcedisplacement transfer function in axial relaxation imaging, including the assumed 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 responses 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, h r (t ), is obtained according 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 re31  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, f P Ri , 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 estimate 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 different 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 linear viscous element, η, as shown in Fig. 2.10. The relationship between stress, σ, and strain, ε, for each of these elements is σe = µε  (2.5)  σv = ηε̇,  (2.6)  and  where the subscripts e and v denote the elastic and viscous components, respectively. 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 relation to the transmission medium. Probe  Transmission Medium  Focal Point Voigt Model  η  µ  33  2.3. Data Processing Software  Figure 2.11: Nyquist plot of normalised Voigt model transfer function, HV (ω). Nyquist plot of normalised Voigt model 0 -0.05 -0.1 -0.15  Imaginary  -0.2 -0.25 -0.3 -0.35 -0.4 -0.45 -0.5  0  0.1  0.2  0.3  0.4  1D Voigt model is HV (ω) =  0.5 Real  0.6  0.7  0.8  0.9  1 ε̂ (ω) = , σ̂ (ω) µ + j ωη  1  (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, f cutoff =  1 µ , 2π η  (2.9)  and that the phase response is −1  ∠ HV (ω) = tan  µ  ¶ −ωη . µ  (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 0 , say. Thus Hr (ω) = K 0  1 , µ + j ωη  (2.11)  where Hr (ω) is the relative transfer function estimated in axial relaxation imaging. 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 =  µ 1 = . 0 K Hr (ω = 0)  (2.12)  The relative viscous parameter, η r , is related to the imaginary component of the relative transfer function: © ª 1 ℑ Hr∗ (ω) η , ηr = 0 = K ω |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 (ω), particularly 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 account 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 and n 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 0 , 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 magnitude 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 parameters µ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 linear 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: Shear Wave 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 using the time of peak displacement in the responses at different lateral locations (the lateral relaxation response). This method or similar has been applied in numerous 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 relaxation 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 location acquired. These estimates are denoted t pk1 , t pk2 , . . . , t pkQ , where Q is the total number of lateral relaxation responses acquired and x 1 , x 2 , . . . , 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 latter 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, lowpass, 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 filter 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 t pk 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 t pk 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, c s . If responses were acquired in only two locations (Q = 2), then the estimate can be simply calculated as cs =  x2 − x1 , t pk2 − t pk1  (2.16)  assuming x 1 is closer to the pushing location than x 2 . 37  2.3. Data Processing Software In estimating c s , 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 s  cs =  E , 2 (1 + ν) ρ  (2.17)  where ν is the Poisson’s ratio of the medium and ρ is the density [36]. For gelatine and soft tissue media, it can generally be assumed ν = 0.5 (implying incompressibility) and ρ = 1000 kg/m3 , giving the equation E = 3000c s2 .  (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 conversion and display occur offline after acquisition of the RF echo signals. However, potential exists for processing and display to occur in near real time, possibly in conjunction with other ultrasound imaging modes, such as B-mode.  38  2.4. Transmission Medium  2.4 Transmission Medium 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. Phantoms 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 gelatinebased phantoms. These phantoms are described in this section. Gelatine has long been used in the construction of tissue-mimicking phantoms 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 tissues. Increasing the concentration of the additive increases the acoustic attenuation (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 glycerol 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 concentrations 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 E g = 0.0034C g2.09 ,  (2.19)  where E g is the Young’s modulus (in kPa) of the gelatine phantoms in that work and C g is the concentration (in g/L). Note that measurement of E g can be vary with boundary conditions and test methods. Nonetheless, (2.19) is used for predicting elastic values of the phantoms in this work. Gelatine phantoms are commonly employed in ARF works. The Duke University group created very soft gelatine phantoms using the method of [20], producing phantoms with values of E g from 0.1–1.6 kPa [40]. This method is also used in [28]. The sonorheometry group from the University of Virginia has developed very soft acrylamide-based phantoms to model the acoustic and mechanical 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 Phantom Construction The phantoms used in this work were made with ~300 bloom, type A, porcineskin gelatine (Sigma-Aldrich Corp., St Louis, MO) and type 50 Sigmacell cellulose (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 continuously. 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 concentrations. A stirring hot plate (PC-420D, Corning Inc., Lowell, MA) was used for heating and stirring, along with a magnetic stirring bar (Fisher Scientific, Ottawa, 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 limited 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 gelatine concentrations produced phantoms that were very soft. This was required to ensure significant displacements could be generated, as the SonixMDP produces 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 plastic 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 oc41  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 University 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 significant effect. As the medium heats during ultrasound imaging, the speed of sound may increase or decrease, causing objects to appear nearer or further, respectively, 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, k c , which can be positive or negative, depending on the medium. I.e. ∆c = k c ∆T,  (2.20)  where ∆c is the change in the speed of sound, k c 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, k c = −2.93 m/s/°C, cow liver, k c = 1.48 m/s/°C, and a rubber phantom, k c = −3.8 m/s/°C [61]. In water and most tissues, k c is positive, though it is negative for fatty tissues [60]. The simulations of [49] gave a maximum displacement error generated by an ARFI imaging sequence of approximately 2.5 µm, using the above k c 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 estimation, displacements in the medium were observed. These displacements occurred 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 appear 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 explanation. 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 k c . 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 displacements 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. Displacement (µm)  Displacement time series at depth of 20 mm  5  60  60  50  50  40  40  Depth (mm)  15 30 20 20 25 10  Displacement (µm)  10  30  20  10  30 0 35 -10 0  2  4  6 Time (s)  8  10  0  -10  0  2  4  6 Time (s)  8  10  12  PVC phantoms, and no ARF displacements were being generated. This was revealed 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 approximately 10 ms, was observed.  2.5.2 Generating Pushing Pulses Another major issue in developing ARF imaging methods concerned the generation of suitably strong pushing pulses. The SonixRP ultrasound platform (Ultrasonix 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. Unfortunately, this switching caused unpredictable corruption of the subsequent RF echo signals used to image the medium. One to four RF signals would be corrupted, either consisting only of zero readings, showing large oscillations like the continuous-wave ARF pushing pulses themselves, or showing large amplitude spikes throughout the signal. An attempt was made to insert “dummy” imaging pulses between the continuous-wave pushing pulse and the first de44  2.5. Development Issues  Figure 2.14: A spatiotemporal displacement profile (left) showing ARFgenerated displacements in a very soft gelatine phantom. An obvious peak occurs at the focal depth of 20 mm, the displacement time series for which is also shown (right). Displacement (µm)  Displacement time series at depth of 20 mm  0  5  5  4.5  5 4  4  10 3.5  Depth (mm)  15 2 20 1  25  Displacement (µm)  3  3 2.5 2 1.5  30  0 1  35 40  -1  0  2  4  6 Time (s)  0.5 0  8 -3  x 10  0  0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 Time (s)  0.01  sired imaging pulse. These pulses were intended to absorb the corrupted signals 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 instead 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 standard pulses (i.e. those generated with trex disabled) was limited to 1.2 µs. Additionally, changing the PRF or pulse duration during a sequence resulted in corruption of RF echoes similar to that which was seen with the trex parameter enabled. Despite numerous experiment with Texo parameter settings, these issues 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 experienced when exceeding a PRF of approximately 20 kHz or increasing the pushing transmit aperture beyond 32 elements. The problem was investigated in experiments 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 ultrasound imaging. An experiment to measure these parameters for the ARF imaging 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 pulseecho ultrasound imaging [64]. These can be separated into thermal and mechanical effects, with the most important effects of each being heating and cavitation. These effects are addressed by different parameters and limits defined in the U.S. Food and Drug Administration (FDA) guidance document for diagnostic 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 occurs 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 probetissue 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 absorption properties of the tissue and the time-averaged intensity of the delivered ultrasound pulses. The spatial-peak, temporal-average intensity parameter, I SP T A , 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 collapse or resonate, releasing significant energy and potentially damaging the tissue in which they formed. The threshold for the occurrence of cavitation is a function of the peak rarefactional pressure, p r , and centre frequency, f c , of the delivered ultrasound pulse [64]. The mechanical index parameter, M I , a unitless value calculated from p r and f c , is given a limit in the FDA guidance document 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, I SP T A and M I , 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 maximum. Each pulse can be described by a pressure-time waveform consisting of a number of sinusoidal cycles of frequency f c , with a total pulse duration of t PD , as shown in the upper part of Fig. 3.1. The time between each pulse is the PRP, TP R , the inverse of the PRF, f P R . The peak rarefactional pressure, p r , is important for determining the mechanical index, M I . The intensity waveform, shown 48  3.1. Background  Figure 3.1: The pressure and intensity waveforms for a fixed-PRF pulse sequence, shown with the important acoustic parameters from the FDA guidance document. Pressure pr tPD TPR  Intensity  PII ISPPA  ISPTA  in the lower part of Fig. 3.1, is derived from the pressure waveform by I=  p2 , ρ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, P I I . The time-average of the intensity over the pulse duration, shown by the upper dotted line, is the spatial-peak, pulse-average intensity, I SP P A . The time-average of the intensity over the PRP, shown by the lower dotted line, is the spatial-peak, time-average intensity, I SP T A . Values of p r , I SP P A and I SP T A are typically measured in a water tank environment. However, delivering equivalent ultrasound pulses into tissue will yield lower values, due to the significantly higher attenuation of tissue, which dissipates 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 applied, 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] ¡ ¢ p r.3 = p r exp −0.0345 f c z SP ,  (3.2)  where z SP is the axial coordinate of the spatial location of the peak rarefactional pressure. Note that in the exponent, f c must be expressed in MHz and z SP in cm. Similarly, applying the derating to the spatial-peak, time-average intensity gives [64] ¡ ¢ I SP T A.3 = I SP T A exp −0.0691 f c z SP ,  (3.3)  where the exponential attenuation factor doubles because intensity is proportional to the square of pressure. The mechanical index can then be calculated as [71] p r.3 MI = p . fc  (3.4)  The intensity parameters are all interrelated [64]: P I I = I SP P A· t P D = I SP T A· TP R .  (3.5)  The limits imposed by the FDA are M I < 1.9, I SP P A.3 < 190 W/cm2 and I SP T A.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 experiment. Parameter Pulse duration, t P D PRF, f P R Transmit voltage, VT X Centre frequency, f c Focal depth, D f Transmit aperture, A T X  Value 0.42 µs 5 kHz 47 V 9.5 MHz 20 mm 32 elements  3.2.1 Method A water tank acoustic intensity measurement system with a membrane hydrophone (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 elements 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 interpolated 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 P I 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 I SP P A and I SP T A 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, P I I , measured in a water tank during the acoustic output experiment. Pulse Intensity Integral  70  0.14  65  0.12  55  0.1  50 45  0.08  40 0.06  35 30  0.04  25 20  0.02  10  8  6  4  2  0  -2  -4  -6  -10  -8  15 10  mJ/cm2  Axial Depth (mm)  60  0  Lateral Location (mm)  experiment, and the experimental value of P I 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 of P I 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 axiallateral scan was performed at this elevation. Measurements of P I I were taken at each scanning step to give a 2D map, shown in Fig. 3.2. The peak in P I 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 P I I , M I and I SP P A . 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 Experiment Pushing Imaging  MI 0.41 0.48 0.48  P I I (mJ/cm2 ) 0.096 1.64 0.034  I SP P A.3 (W/cm2 ) 61.2 90.1 90.1  I SP T A.3 (mW/cm2 ) 129 13000 135  values and the interpolated values for typical ARF pushing and imaging pulses are given in Table 3.2. Note that the value of M I 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 M I are well below the limit of 1.9. Similarly, the values of I SP P A.3 are well below the 190 W/cm2 limit. These results imply that cavitation will not occur in tissue undergoing ARF imaging with the system of this work, as the instantaneous pressures generated are not high enough. Clearly, the interpolated value of I SP T A.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 I SP T A.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 sequence 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 I SP T A.3 of a similar order to the results here [40, 42]. For example, they measure 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 I SP T A.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 simulations confirm this is the case for most situations [49]. Works regarding the SSI technique [2, 65], approach the issue of the applicability of the I SP T A.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 intensity imaging pulses. In assessing the safety of the method, they effectively calculate the P I 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 I SP T A.3 is calculated. Applying this concept to typical sequences in this work, which use 100 pushing pulses and 300 imaging pulses per sequence loop, the values of P I 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 I SP T A.3 . Although this argument may be reasonable, 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 overheating 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 techniques 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 platform to the probe. The energy delivered to each element depends on the pulse transmit voltage, 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 consists 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, TP Ri . The pushing phase is a series of N p pushing pulses (labelled p.1 to p.N p ) transmitted at the pushing PRP, TP Rp . Each imaging pulse consists of M i cycles at the centre frequency, f c , while each pushing pulse consists of M p 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, E T , over L repetitions of the sequence, assuming a 50 W probe impedance, is given by ET =  VT2 X M i Ni + M p N p 100  fc  L.  (3.6)  This quantity represents the power delivered by the transmitted signal (VT2 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. Mi cycles (tPDi)  Mp cycles (tPDp)  Imaging pulse  VTX  Pushing pulse  T c =1/f c  i.1  i.2  T PRi =1/f PRi Imaging phase  . . .  i.Ni  p.1  p.2  . . .  p.Np  i.1  T PRp =1/f PRp Pushing phase  Sequence loop  is looped, t T . That is, ¡ ¢ t T = Ni TP Ri + N p TP Rp 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 ultrasound 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 ability 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 sequence 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, E T,M AX , is allowable be-  56  3.3. Probe Heating and Protection fore damage to the transducer occurs, then the maximum number of sequence repetitions, L M AX , can be determined from (3.6). I.e. L M AX = E T,M AX  100 VT2 X  fc . M i Ni + M p N p  (3.8)  A suitable value for E T,M AX 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, f c = 6.7 MHz, M i = 1, Ni = 300, M p = 48, and N p = 100, giving a total energy, E T = 0.589 J. Considering VT X will remain at 47 V in all ARF experiments of this work, (3.8) can be re-written in a form suitable for programming into the ARF Imaging program: L M AX <  fc ¡  38 M i Ni + M p N p  ¢,  (3.9)  where f c is in Hz. In the Texo program, the value of L M AX 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 L M AX 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 suitable, 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 E T,M AX would need to be determined.  57  Chapter 4  ARF Experiments 4.1 Introduction This chapter details experiments performed to confirm some of the theory related to ARF imaging, verify that the ultrasound system was capable of generating ARF-induced displacements in soft tissue-like transmission media, and explore 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 Section 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 measure 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 parameters examined were: the number of cycles per pushing pulse, M p ; the pushing PRF, f P R p ; the number of pushing pulses per sequence loop, N p ; the centre frequency, f c ; and the pushing transmit aperture, A T X p . Table 4.1 outlines the values that were tested for each parameter. Note that the values selected for f c 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 Mp fP R p Np fc AT X p  Values 12, 24, 36, 48 cycles 10, 20, 30,..., 80 kHz 25, 50, 75, 100, 125, 150 pulses 4.0, 5.0, 6.7, 8.0, 10.0, 13.3 MHz 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 parameter 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 (M p = 48, f P R p = 20 MHz, N p = 100, f c = 6.7 MHz and A T X p = 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 – M p , f P Rp , N p , f c and A T X p – are shown in Figs. 4.3–4.7. In all results, excepting those where significant displacement was not observed (e.g. the f c 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 f P Rp = 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 response at the focal depth (right). Displacement (µm)  Relaxation response at focal depth 30 25  10  25  20  20 15 25 10  Displacement (µm)  Depth (mm)  15  20  15  10  30 5 5  35  0  1  2 3 Time (ms)  4  0  0  2  4  6  8  10 12 Time (ms)  14  16  18  20  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 A T X p 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 parameters, 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 exponential 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 M p , the number of cycles per pushing pulse, show that peak displacement, and therefore ARF, increases with M p . The trend is ini61  4.4. Discussion  Figure 4.3: The results for the experiment in which the number of cycles per pushing pulse, M p 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. Relaxation responses at focal depth for various values of Mp  Peak displacements for various values of Mp  30  30 12 cycles 24 cycles 36 cycles 48 cycles  25 Peak displacement (µm)  Displacement (µm)  25  20  15  10  15  10  5  0  20  0  0.5  1  1.5  2  2.5 3 Time (ms)  3.5  4  4.5  5 10  5  15  20 25 30 35 40 Number of cycles per pushing pulse, Mp  45  50  Figure 4.4: The results for the experiment in which the pushing PRF, f P Rp , 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. Relaxation responses at focal depth for various values of pushing PRF  Peak displacements for various values of pushing PRF 38  40 10 kHz 20 kHz 30 kHz 40 kHz 50 kHz 60 kHz 80 kHz  Displacement (µm)  30 25  36 34 Peak displacement (µm)  35  20 15  32 30 28 26 24  10 22 5 0  20  0  0.5  1  1.5  2  2.5 3 Time (ms)  3.5  4  4.5  5  18 10  20  30  40 50 Pushing PRF (kHz)  60  70  80  62  4.4. Discussion  Figure 4.5: The results for the experiment in which the number of pushing pulses per sequence loop, N p , 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. Relaxation responses at focal depth for various values of Np  Peak displacements for various values of Np  30  29 25 pulses 50 pulses 75 pulses 100 pulses 125 pulses 150 pulses  20  28 27 Peak displacement (µm)  Displacement (µm)  25  15  10  26 25 24 23 22  5 21 0  0  0.5  1  1.5  2  2.5 3 Time (ms)  3.5  4  4.5  20 20  5  40  60 80 100 120 Pushing pulses per sequence loop, Np  140  160  Figure 4.6: The results for the experiment in which the transmit centre frequency, f c , 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. Peak displacements for various values of centre frequency  Relaxation responses at focal depth for various values of centre frequency 35 4 MHz 5 MHz 6.7 MHz 8 MHz 10 MHz 13.3 MHz  Displacement (µm)  25  30  25 Peak displacement (µm)  30  20 15 10  20  15  10  5 5 0 -5  0 0  0.5  1  1.5  2  2.5 3 Time (ms)  3.5  4  4.5  5  4  5  6  7 8 9 10 11 Centre frequency, fc (MHz)  12  13  63  4.4. Discussion  Figure 4.7: The results for the experiment in which the pushing transmit aperture, A T X p , 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. Peak displacements for various values of pushing transmit aperture 30 28 26 Peak displacement (µm)  Displacement (µm)  Relaxation responses at focal depth for various values of pushing trans mit aperture 30 16 elements 32 elements 48 elements 25 64 elements 20  15  10  24 22 20 18 16  5 14 0  0  0.5  1  1.5  2  2.5 3 Time (ms)  3.5  4  4.5  5  12 15  20  25  30 35 40 45 50 55 Pushing transmit aperture (no. elements)  60  65  tially linear, with the displacement doubling from 9 to 18 µm as M p doubles from 12 to 24. This agrees with the basic theory of ARF, particularly (1.1). The timeaverage acoustic intensity in this equation, 〈I 〉, is proportional to M p (M p essentially 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 M p . After these initial results, the peak displacement increases more slowly as M p 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 M p . This effect seems to be apparent in the other results of this chapter also. In any case, the maximum value, M p = 48 cycles, which is the maximum allowed by the SonixMDP firmware, gives the largest displacement, and therefore is the most suitable for ARF generation. In a similar result is observed in Fig. 4.4 for pushing PRF, f P Rp . Peak displacement generally increases with f P Rp , 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, A T X p . Displacement (µm) for Ap = 32 elements  Displacement (µm) for Ap = 16 elements 15 10  10 20 25 5  Depth (mm)  Depth (mm)  15  10  25  15  20  20  15  25  10  30  30  5 35  35 0 0  5  10 Time (ms)  0 0  15  Displacement (µm) for Ap = 48 elements  5  10 Time (ms)  15  Displacement (µm) for Ap = 64 elements 14  10  20 15  20 10  25 30  5  35  Depth (mm)  Depth (mm)  15  10  12  15  10 8  20  6 25 4 30  2  35  0  0 0  5  10 Time (ms)  15  0  5  10 Time (ms)  15  65  4.4. Discussion at higher values of f P R p . The peak displacement eventually saturates at 60 kHz, and even reduces at 80 kHz. Following the argument for M p above, time-average intensity, ARF and peak displacement should all be proportional to f P Rp , because f P R p also controls the pushing pulse duty cycle. Again, the high levels of power drawn from the SonixMDP power supply at high values of f P Rp 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 M p and f P R p , the number of pushing pulses per sequence loop, N p , does not affect the time-average intensity and the magnitude of the ARF generated. Instead, it affects the total time for which the ARF is applied. This allows 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 N p in Fig. 4.5 shows that the peak displacement gradually saturates over a range from 25–100 pulses per sequence loop, matching 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 N p increase the peak displacement, eventually increasing N p will have little effect. The reduction in peak displacement seen after N p = 100 in Fig. 4.5 is again likely due to the power supply issues in the SonixMDP. Thus, N p = 100 is an optimal default value, as it ensures saturation of displacement without over-straining the power supply. Displacement 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, f c , experiments essentially reveal the bandwidth of the L14–5/38 ultrasound probe. This probe has a nominal bandwidth 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 frequencies, 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, A T X p , 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 additional power draw of the added elements. Thus, A T X p = 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 ultrasound, 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 regarding 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 homogeneous 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 imaging 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 ultrasound 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 representing 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 phantoms. The values of K and K 0 (see Equations (2.1) and (2.11)) could therefore be assumed to be equal for each phantom, as is necessary for comparing relative 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 water 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 relaxation imaging, relaxation responses were acquired at six lateral locations spaced every 0.9144 mm from x 1 = 1.829 mm to x 6 = 6.401 mm. This corresponds to a shift of three probe elements between each location. Note that x = 0 corresponds 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 phantoms 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) centred 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. Displacement (µm) for phantom A1  Displacement (µm) for phantom A2  10  60  10 50 45  15  50  15  40  30 25  25  20  20 Depth (mm)  Depth (mm)  40  35  20  30 25 20  15  30  30  10  10  5  35  35  0 0  5  10  15 Time (ms)  20  0  25  0  5  Displacement (µm) for phantom B1  10  15 Time (ms)  20  25  Displacement (µm) for phantom B2  10  10 20  20  15  15 15  10  25  30  Depth (mm)  Depth (mm)  20  20  15  25  10  30  5  5  35  35 0 0  5  10  15 Time (ms)  20  0  25  0  5  Displacement (µm) for phantom C1  10  15 Time (ms)  20  25  Displacement (µm) for phantom C2 12  10  12  10  10  10  15  15 8  8 20  6 25 4  Depth (mm)  Depth (mm)  20  6 25  30  2  30  35  0  35  4  2  0  -2  -2 0  5  10  15 Time (ms)  20  25  0  5  10  15 Time (ms)  20  25  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. Focal-depth displacement relaxation responses 60 Phantom A Phantom B Phantom C  50  Displacement (µm)  40  30  20  10  0  -10  0  5  10  15 Time (ms)  20  25  30  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 relaxation imaging transfer function estimation results of Section 5.3.2 gave estimates 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. Nyquist plots of relative force-displacement transfer function estimates 0.1 1000 Hz 0  0 Hz  0 Hz  0 Hz  -0.1  Imaginary  100 Hz -0.2  -0.3  300 Hz  -0.4 100 Hz Phantom A Phantom B Phantom C  -0.5  -0.6 -0.2  0  0.2  0.4  0.6  0.8  1  1.2  Real  Figure 5.4: The magnitude (left) and phase (right) plots of the relative forcedisplacement transfer function estimates for the three phantom types. Results are normalised to the mean of the 0 Hz values of the type-A phantoms. Magnitude plots of relative force-displacement transfer function estimates  Phase plots of relative force-displacement transfer function estimates 0  0  Phantom A Phantom B Phantom C  Phantom A Phantom B Phantom C  -0.5  -1 Phase (rad)  Relative Magnitude  10  -1  10  -1.5  -2  -2.5  -3  2  2  10  10  Frequency (Hz)  Frequency (Hz)  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 typeA 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). Relative viscosity estimates with frequency  Relative elasticity parameter estimates 4  5  Phantom A  4.5  Phantom B Phantom C  3.5  4  Pow er Law  3 Relative Viscosity  Relative Elasticity  3.5 3 2.5 2  2.5 2 1.5  1.5 1 1 0.5  0.5 0  A1  A2  B1  B2 Phantom  C1  C2  0  0  100  200  300 400 500 Frequency (Hz)  600  700  800  τ. 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 results, 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 predicted 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 using a linear least-squares fit to logarithmic plots of the relative viscous parameter 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 n A = −0.10, n B = −0.14 and nC = −0.18 for phantom types A, B and C, respectively. The relaxation 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. Time constant estimates with frequency 2.4 Phantom A  2.2  Phantom B Phantom C  2  Time Constant (ms)  1.8 1.6 1.4 1.2 1 0.8 0.6 0.4  0  100  200  300 400 500 Frequency (Hz)  600  700  800  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 relaxation 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 resulting 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). Lateral relaxation responses at focal depth for type-A phantoms 30 1.88 mm 2.82 mm 3.77 mm 4.71 mm 5.65 mm 6.59 mm  25  Displacement (µm)  20  15  10  5  0  -5  0  5  10  15 Time (ms)  20  25  30  Lateral relaxation responses at focal depth for type-B phantoms 12 1.88 mm 2.82 mm 3.77 mm 4.71 mm 5.65 mm 6.59 mm  10  Displacement (µm)  8  6  4  2  0  -2 0  5  10  15 Time (ms)  20  25  30  Lateral relaxation responses at focal depth for type-C phantoms 8 1.88 mm 2.82 mm 3.77 mm 4.71 mm 5.65 mm 6.59 mm  7 6  Displacement (µm)  5 4 3 2 1 0 -1 -2  0  5  10  15 Time (ms)  20  25  30  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. -3  9  A1 A2 B1 B2 C1 C2  8 7 Time of peak displacement (s)  Lateral location vs time of peak displacement  x 10  6 5 4 3 2 1 0  1  2  3  4 5 Lateral Location (m)  6  7 -3  x 10  5.4.2 Shear Wave 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, c s , 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, E r , were also calculated by dividing by the mean absolute value of the type-A phantoms. Absolute and relative Young’s modulus values predicted by the empirical formula (2.19), E g and E g ,r , respectively, were calculated for comparison. All these results 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 (c s ) and Young’s modulus (E ) estimation results of lateral relaxation imaging on the six homogeneous phantoms. The empirically predicted Young’s moduli are also shown (E g ). Relative Young’s modulus or elastic parameter values from lateral relaxation imaging (E r ), the empirical formula (E g ,r ) and axial relaxation imaging (µr ) are given in the last three columns for comparison. For each estimation method, relative values are normalised by the mean value of the type-A phantoms. Phantom A1 A2 B1 B2 C1 C2  c s (m/s) 0.70 0.69 1.04 1.06 1.48 1.41  E (kPa) 1.46 1.42 3.24 3.39 6.61 6.00  E g (kPa) 1.78 1.78 4.16 4.16 7.58 7.58  Er 1.01 0.99 2.25 2.36 4.60 4.17  E g ,r 1 1 2.34 2.34 4.26 4.26  µr 1.03 0.97 2.35 2.26 4.63 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 behaviour. 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 acquisitions, 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 de77  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, indicating 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 magnitude and phase plots, a low-pass response is observed, agreeing with the Voigt model. Characteristics that distinguish between phantoms as gelatine concentrations decrease include a larger 0 Hz magnitude component, a steeper magnitude 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 parameter 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 parameter, as mentioned in Section 2.4.1. These results show that parameters derived directly from the transfer function estimates could be useful in viscoelastic characterisation 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 values 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 expected, 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 decrease 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 models 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 constant 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 values 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 successfully in characterisation of soft tissue.  5.5.2 Lateral Relaxation Imaging All the lateral relaxation responses of Fig. 5.7 clearly reveal the lateral propagation 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 relaxation responses generated with the 100 Hz component of the relaxation time constant estimate. Generated responses are assumed to have an initial displacement equal to that of the observed responses. Experimental and time constant estimate relaxation res ponses 60 Phantom A Phantom B Phantom C Time Constant Results  50  Displacement (µm)  40  30  20  10  0  -10  0  5  10  15 Time (ms)  20  25  30  location have a peak displacement that occurs later in time and with a lower amplitude, 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 ab80  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 maximum error of 10.8% between any two methods. Finally, by assuming the elasticity parameter of the Voigt model is approximately 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, respectively. 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 instead 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 imaging that gives the ability to extract frequency-dependent viscoelastic information 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 relaxation 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 repetitions 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 limiting 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 system 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 observed to create a near-saturated displacement, validating the assumption in axial 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 empirical predictions and lateral relaxation imaging estimates of relative Young’s modulus. It was also shown that the relaxation time constant parameter gave good separation between phantom types, without being dependent on the magnitude 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 viscoelastic parameter estimation: A novel new method of estimating the forcedisplacement 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 parameter is further modelled with a power law, giving a flow index parameter estimate. Considering the range of information obtained, the acquisition 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 ultrasound 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 single 5 ms pushing phase. Periodic, swept-frequency methods use continuous wave ARF, which may be applied for over 500 ms to obtain an equivalent frequency response (e.g. [28]). The shorter excitation and monitoring 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 researchers to follow this work and develop their own ARF imaging techniques. 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 overhead of work in the construction of a custom ultrasound system. The description 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 Future Work The ARF imaging techniques presented in this thesis are relatively new developments, 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 currently 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 generate detectable displacements in stiffer media. Again, this may be achieved with external hardware. Note that it would require re-evaluation of acoustic 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 relation 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 relaxation imaging technique simply mimics other techniques from the literature. Initial work has begun on developing transfer function estimation 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 characterisation. • 2D axial relaxation imaging: More work is required to develop 2D axial relaxation imaging methods for heterogeneous media. The spatial variability 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 imaging of soft tissue shear viscosity,” in IEEE Ultrasonics Symposium, Honolulu, HI, October 2003, pp. 937–940. [2] J. Bercoff, T. Mickaël, and M. Fink, “Supersonic shear imaging: A new technique 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 Control, 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 Transform and Its Applications, 3rd ed. New York: McGraw-Hill, 2000, ch. 4, p. 63. [6] Canadian Cancer Society’s Steering Committee, “Canadian cancer statistics 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, December 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 measuring tissue elasticity and viscosity,” IEEE Transactions on Ultrasonics, Ferroelectrics 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,” IEEE Transactions 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 estimation based on spectral analysis,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 55, no. 7, pp. 1611–1625, July 2008. [15] M. Fatemi and J. F. Greenleaf, “Ultrasound-stimulated vibro-acoustic spectrography,” 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 in Medicine and Biology, 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 characterization 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 Frequency Control, vol. 44, no. 6, pp. 1355–1365, November 1997. [21] C. Joly-Duhamel, D. Hellio, and M. Djabourov, “All gelatin networks: 1. Biodiversity 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 in Medicine 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 Biology, vol. 29, no. 10, pp. 1405–1413, October 2003. [24] E. E. Konofagou, J. D’hooge, and J. Ophir, “Myocardial elastography—a feasibility 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, “Elastic 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 resonance 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 Acoustical 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 Biomedical Imaging, Arlington, VA, April 2007, pp. 892–895. [29] D. Liu and E. S. Ebbini, “Viscoelastic property measurement in thin tissue constructs using ultrasound,” IEEE Transactions on Ultrasonics, Ferroelectrics 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 mimicking 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 heterogeneous 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 acoustical imaging techniques for nondestructive evaluation of art objects,” Research 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 lesions: Preliminary results in the breast,” Ultrasound in Medicine and Biology, 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 response 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 feasibility,” Ultrasound in Medicine and Biology, vol. 28, no. 2, pp. 227–235, February 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 excitation in liver in vivo,” in IEEE Ultrasonics Symposium, Vancouver, BC, October 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,” Ultrasonic 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 Biology, 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 impulsive acoustic radiation force,” IEEE Transactions on Ultrasonics, Ferroelectrics 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. Nightingale, “Dynamic mechanical response of elastic spherical inclusions to impulsive acoustic radiation force excitation,” IEEE Transactions on Ultrasonics, 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 microscopic 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 vessel 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 response 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 duality 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. Montaldo, M. Muller, A. Tardivon, and M. Fink, “Quantitative assessment of breast lesion viscoelasticity: Initial clinical results using supersonic shear imaging,” Ultrasound in Medicine and Biology, vol. 34, no. 9, pp. 1373–1386, September 2008. [66] Ultrasonix Medical Corporation, “Transducer specification sheet,” Ultrasonix Medical Corp., August 2009. [67] Ultrasonix Medical Corporation, “Texo documentation 2.0,” Ultrasonix Medical Corp., July 2010. [Online]. Available:  http://ultrasonix.com/  wikisonix/api/texo/2.0/ [68] M. W. Urban, S. Chen, and J. F. Greenleaf, “Error in estimates of tissue material 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 Transducers, U.S. Food and Drug Administration Std., September 2008. [Online]. Available:  http://www.fda.gov/downloads/MedicalDevices/  DeviceRegulationandGuidance/GuidanceDocuments/UCM070911.pdf [72] J. Vappou, C. Maleke, and E. E. Konofagou, “Quantitative viscoelastic parameters measured by harmonic motion imaging,” Physics in Medicine 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 potential,” 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 radiation 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. Hobson, “Investigating the effects of viscosity on focused, impulsive, acoustic radiation force induced shear wave morphology,” in IEEE Ultrasonics Symposium, 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 on Ultrasonics, Ferroelectrics 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 radiation 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 volume, 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 general 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 assumed form of the wave and medium, is given by ¶ ∂v ∂v ∂p 4 ∂2 v ρ +v =− + (η0 + η) 2 . ∂t ∂z ∂z 3 ∂z µ  (A.1)  Here, ρ is the density of the medium, v is velocity, t is time, z is the spatial coordinate for the direction of propagation (assumed to be axial), p is pressure, and η0 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 4 ∂2 v + (η0 + η) 2 . ∂z 3 ∂z  (A.2)  The equation of continuity, ¡ ¢ ∂ρ ∂ ρv + = 0, ∂t ∂z  (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 ∂p 4 ∂2 v + ρv +v =− + (η0 + η) 2 . ∂z ∂z ∂z 3 ∂z  (A.4)  Here, the last two terms of the left-hand side represent additional inertial components resulting from the fact that the medium has not been assumed to be incompressible. Together, they constitute an excess force per unit volume term, −F E , such that  ¡ ¢ ∂ ρv ∂v − F E = ρv +v . ∂z ∂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 successively higher-order approximations. Thus, p 0 + ²p 1 + ²2 p 2 + . . .  p  =  ρ  = ρ 0 + ²ρ 1 + ²2 ρ 2 + . . .  v  = ²v 1 + ²2 v 2 + . . . ,  (A.7)  where the subscripts denote the order of the approximations. The zeroth-order values, p 0 and ρ 0 , are the static values of pressure and density (while v 0 = 0). Using zeroth-order (²0 ) terms in (A.4) gives  ∂p 0 ∂z  = 0, such that p 0 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  ∂v 1 ∂p 1 4 ∂2 v 1 =− + (η0 + η) 2 . ∂t ∂z 3 ∂z  (A.8)  This is equivalent to the equation for an incompressible medium (the −F E terms disappear), and describes the linear wave propagation at the fundamental harmonic frequency. The velocity here is assumed to have the form of a damped harmonic travelling wave: v 1 = v p e −αz sin (ωt − kz) ,  (A.9)  where v p 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 ω2 0 α= η + η , 3 2ρ 0 c 3  where c =  q  p0 ρ0  (A.10)  is the speed of sound.  Now considering only the second-order (²2 ) terms, (A.4) becomes, after simplification, ¢ ∂ ¡ 4 ∂2 v 2 ∂v 1 −∂p 2 ρ 0 v 2 + ρ 1 v 1 + 2ρ 0 v 1 = + (η0 + η) 2 . ∂t ∂z ∂z 3 ∂z  (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]):  ¿  À ¢ ∂ ¡ ρ 0 v 2 + ρ 1 v 1 = 0. ∂t  (A.12)  Hence, (A.11) becomes ¿ À ¿ À ∂v 1 −∂p 2 4 ∂2 v 2 0 2ρ 0 v 1 = + (η + η) 2 , ∂z ∂z 3 ∂z  or, ¿  〈−F E 2 〉 =  −∂p 2 4 ∂2 v 2 + (η0 + η) 2 ∂z 3 ∂z  (A.13)  À  ,  (A.14)  99  Appendix A. Derivation of an Expression for the Acoustic Radiation Force where 〈−F E 2 〉 is the time-average of the second-order terms of the excess force per unit volume, F E . The time-average excess force per unit volume is thus ¿ À ∂v 1 〈F E 2 〉 = F = −2ρ 0 v 1 , ∂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ρ 0 v p2 e −2αz α sin2 (ωt − kz) + k sin (ωt − kz) cos (ωt − kz)  = αρ 0 v p2 e −2αz ,  (A.16)  because the time-average of the sine-squared term is α/2, while that of the sinecosine 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 = ρ 0 cv [10]. Thus, 〈I 〉 = =  ­  ρ 0 c v 12  ®  ρ 0 c v p2 e −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 constant, 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 )]  ˆ∞  h (t ) d t − K θ (t ) ⊗ h (t )  = K −∞  = K H (ω = 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) = K H (ω = 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 0 d t 0,  θ (t ) ⊗ f (t ) =  (B.3)  −∞  for a function f (t ). Using (B.1) and (B.3), gives  ˆ0 ¡ ¢ h t0 dt0  K s (0) = K −∞  = 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) = K H (ω = 0). Then, from (B.2), r (t ) = r (0) − K s (t ) and K s (t ) = r (0) − r (t ), as required.  102  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items