A Wave Equation Approach to Ultrasound Elastography by Ali Baghani B.A.Sc., Sharif University of Technology, 2004 M.A.Sc., University of Tehran, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University of British Columbia (Vancouver) January, 2010 c Ali Baghani 2010 Abstract Tissue elastography is a field of medical imaging which deals with obtaining images of tissue mechanical properties, in particular tissue elasticity. Every elastography imaging system has three major components: an exciter, an imaging system, and an image processing software module. The principle of operation of these systems is that softer and harder tissues react differently to a mechanical excitation. The imaging system captures images of the tissue as the mechanical exciter deforms the tissue. The image processing software module then computes the tissue motion from the acquired images, and solves an inverse problem to obtain tissue elasticity from tissue displacements. In ultrasound elastography all three components of the displacement are not available, and the measurement is carried out only on the imaging plane. The problem of inversion using partial measurements is an inherently challenging problem. To solve this problem the phase speed is estimated in the first place. This thesis studies the process of estimating tissue elasticity from tissue displacements, using the phase speed estimation. It is shown that, in general, the elasticity cannot be inferred from the phase speed, but under certain boundary and excitation conditions, a relationship exists between the phase speed and the elasticity. However the relationship depends on the boundary and excitation conditions. Based on the developed results, a novel rheometry method is proposed for visco-elastic characterization of tissue samples. To increase the resolution of the elasticity images it is of interest to use higher frequencies of excitation. However the inherent low frame rate of ultrasound systems posed a limitation. A novel high-framerate ultrasound system is introduced which is capable of tracking motions of up to 500Hz. The high-frame-rate system is used in conjunction with the inversion algorithms to form an elastography system. The performance of the system is tested experimentally on tissue mimicking material having hard and soft inclusions. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Abbreviations Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . xii . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . 1.1 Motivation . . . . 1.2 Background . . . 1.3 Thesis Objectives 1.4 Chapter Summary References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 5 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 Measurement of Viscoelastic Properties of Tissue Mimicking Material Using Longitudinal Wave Excitation . . . . . 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Theory of Wave Motion . . . . . . . . . . . . . . . . . . . . . 2.2.1 Longitudinal Plane Waves in an Infinite Medium . . . 2.2.2 Longitudinal Plane Waves in a Finite Medium . . . . 2.3 Experimental Methods and Setup . . . . . . . . . . . . . . . 2.4 Proposed Rheometry Methods . . . . . . . . . . . . . . . . . 2.4.1 Peak to Node Method . . . . . . . . . . . . . . . . . . 2.4.2 Model Fitting Method . . . . . . . . . . . . . . . . . 2.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 15 17 19 20 22 25 28 29 30 iii Table of Contents 2.6 2.7 2.5.1 Peak to Node Results 2.5.2 Model Fitting Results Discussion . . . . . . . . . . Summary and Conclusion . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 30 33 39 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3 Theoretical Limitations of the Elastic Wave Equation Inversion for Tissue Elastography . . . . . . . . . . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Local Inversion of the Wave Equation . . . . . . . . . . . . . 3.3 A Critique on the Terminology Associated with the Wave Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Phase Speed of Mechanical Waves in Elastic Mediums . . . . 3.4.1 Dependence of the Phase Speed on the Geometry of the Wave: Infinite Mediums . . . . . . . . . . . . . . 3.4.2 Dependence of the Phase Speed on the Geometry of the Medium: Wave Guides . . . . . . . . . . . . . . . 3.5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Shear Waves in a Tissue Mimicking Infinite Medium . 3.5.2 Waves in a Tissue Mimicking Cylindrical Rod . . . . 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 The Influence of the Boundary Conditions on Longitudinal Wave Propagation in a Viscoelastic Medium . . . . . . . . . . . . . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Governing Equations in a Linear Viscoelastic Medium . . . . 4.3 Wave Equation and the Propagation Speed . . . . . . . . . . 4.3.1 Longitudinal Wave in a Finite Medium Approximating a Thin Rod . . . . . . . . . . . . . . . . . . . . . 4.3.2 Longitudinal Wave Induced by a Point Source Exciter 4.4 Determining the Wave Speed from the Spectra of the Displacements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Case (i) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Case (ii) . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Dynamic Finite Element Analysis . . . . . . . . . . . 46 46 49 52 55 55 59 61 61 64 67 70 78 78 80 82 83 84 87 89 90 91 91 iv Table of Contents 4.6 4.7 4.5.2 Simulation of the 4.5.3 Experiments . . Discussion . . . . . . . Conclusion . . . . . . . References Wave Propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5 A High Frame Rate Ultrasound System for the Study Tissue Motions . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Sector-Based High Frame Rate Acquisition Technique . . . 5.2.1 Sequencing . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Time Delay Compensation . . . . . . . . . . . . . . 5.2.3 Implementation Issues . . . . . . . . . . . . . . . . . 5.3 Requirements on the Tissue Motion . . . . . . . . . . . . . 5.4 Applications and Experimental Results . . . . . . . . . . . 5.5 Summary and Conclusion . . . . . . . . . . . . . . . . . . . References of . 107 . 107 . 109 . 109 . 113 . 116 . 119 . 122 . 128 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6 Conclusion . . . . . . . . . . . . . . . . . . . . 6.1 Summary . . . . . . . . . . . . . . . . . . . 6.2 Contributions . . . . . . . . . . . . . . . . 6.2.1 Rheometry . . . . . . . . . . . . . . 6.2.2 Poly Vinyl Chloride Phantoms . . . 6.2.3 Phase Speed and Tissue Elasticity . 6.2.4 High Frame Rate Ultrasound System 6.2.5 Ultrasound Elastography System . . 6.3 Future Work . . . . . . . . . . . . . . . . . References . 93 . 96 . 98 . 101 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 134 136 136 136 136 138 138 139 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 Appendices A Conventional Rheometry . . . . . . . . . . . . . . . . . . . . . 143 B Tissue Mimicking Materials for Ultrasound Elastography 148 B.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . 148 B.2 PVC Phantom Recipe . . . . . . . . . . . . . . . . . . . . . . 152 B.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 152 v Table of Contents B.2.2 Preparation Steps . . . . . . . . . . . . . . . . . . . . 152 B.2.3 Other Important Notes and Tips . . . . . . . . . . . . 153 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 C Phasor Calculations . . . . . . . . . . . . . . . . . . . . . . . . 158 vi List of Tables 2.1 2.2 2.3 2.4 24 24 25 2.6 Fabrication and measured properties of the five test phantoms Parameters of the ultrasound imaging system . . . . . . . . . Parameters of the speckle tracking motion estimation algorithm Estimated values of longitudinal wave speed cfin and Young’s modulus for the test phantoms using the peak to node method Estimated values of longitudinal wave speed cfin and Young’s modulus for the test phantoms using the model fitting method Summary of the elastic properties of the five test phantoms . 3.1 3.2 Mechanical property values used for the simulation . . . . . . Mechanical property values used for the simulation . . . . . . 61 64 4.1 The theoretical and estimated wave speed for three different conditions. The estimated wave speed are based on three method: first, the frequency that the peak of the spectrum happens was used to reconstruct the wave speed; second, the previous value was enhanced to account for the effect of viscosity; third, the wave-front due to a transient step was analyzed to recover the speed. . . . . . . . . . . . . . . . . . . . . The theoretical and estimated wave speed for an experiment with a soft PVC phantom. A flat and a point-source excitation were applied to explore the effect of boundary conditions on the speed of longitudinal waves. . . . . . . . . . . . . . . . 2.5 4.2 5.1 32 33 37 97 98 Parameters of the high frame rate system . . . . . . . . . . . 124 vii List of Figures 2.1 2.2 2.3 2.4 2.5 2.6 2.7 Simple schematic of the rheometry device . . . . . . . . . . . A typical normalized displacement phasor profile, belonging to phantom No.3 at a frequency of 90Hz. The (a) magnitude and (b) phase of the phasor profile (c) the same normalized phasor profile in the complex plane. The profile is normalized to the phasor of the point at x = 40mm. . . . . . . . . . . . . Fitted model to the measured displacement data. The data is the same as that of Fig. 2.2, which is from phantom No. 3 at 90Hz. (a) absolute value and (b) phase of the phasor profile (c) the same normalized phasor profile in the complex plane. Estimated cfin (ω) for phantom No.3. Plotted in dark solid line is the selected region for estimation of average speed. At lower frequencies no peak occurs inside the block length. . . . (a) Estimated value of the Young’s modulus and (b) relaxation times for the five phantoms at different frequencies of excitation. The lower frequency graphs (4–40Hz) are from conventional rheometry and the higher frequency graphs (above 40Hz) are from model fitting. . . . . . . . . . . . . . . . . . . (a) The result of the excitation with a plate which covers the top surface of the phantom (b) The result of the excitation with a point source. . . . . . . . . . . . . . . . . . . . . . . . Obtained Young’s modulus as a function of the volume portion of Plastic Softener or Hardener added during the fabrication of the phantoms . . . . . . . . . . . . . . . . . . . . . . 23 26 31 32 34 36 37 viii List of Figures 3.1 3.2 3.3 3.4 4.1 4.2 4.3 4.4 Infinite medium; shear beam patterns for three different frequencies all sharing the same phase speed of 5m/s. (a), (b), (c) radial components of the displacement (d), (e), (f) zcomponent of the displacement at 40Hz, 65Hz, and 100Hz respectively. Note that only a cubic portion of the medium near the axis of symmetry is shown; The medium extends in all directions to infinity. . . . . . . . . . . . . . . . . . . . . . Infinite medium; shear beam patterns for three different frequencies all sharing the same phase speed of 12m/s. (a), (b), (c) radial components of the displacement (d), (e), (f) z-component of the displacement at 40Hz, 65Hz, and 100Hz respectively. Note that only a cubic portion of the medium near the axis of symmetry is shown; The medium extends in all directions to infinity. . . . . . . . . . . . . . . . . . . . . . The first four modes of the cylindrical bar . . . . . . . . . . . Permissible beam patterns (modes) in an infinite cylindrical rod for waves having a frequency of 100Hz. Three (and only three) phase speeds are possible. (a), (b), (c) radial components of the displacement (d), (e), (f) z-component of the displacement for phase speeds of 3.0339m/s, 3.719m/s, and 5.630m/s respectively. . . . . . . . . . . . . . . . . . . . . . . Displacement vector in the cylindrical coordinate system. . . The magnitude of the displacement in a 1D model of the wave equation with the purely elastic wave speed of 3.16m/s, length of 6cm, and Young’s modulus of 10kP a. (a) and (b) are the magnitude and contour plots for the case in which the viscosity is 5P as. The frequencies denoted by case (i) and case (ii) are almost equal. (c) and (d) show the the plots for the case that the viscosity of the medium is 30P as. In this case, f1 over-estimates the actual value, while f2 results in a relatively small under-estimation. . . . . . . . . . . . . . The geometry of the 3D region used for FEM simulations. . . The magnitude (a) and the phase (b) of the transfer function between the center point of the phantom and the middle point at the top. The solid line shows the case where the four vertical sides of the phantom are free to move and the dashed lines are for the case where the sides of the phantom are confined from bulging. The gray dash-dot line shows the transfer function for a one-node excitation at the middle point. 62 63 64 66 85 90 92 94 ix List of Figures 4.5 4.6 4.7 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 Effect of the viscosity on the magnitude of the transfer function between the center point and the middle point. The excitation was applied to the entire top surface while the sides were free to bulge. . . . . . . . . . . . . . . . . . . . . . . . . 95 Transient step response for different boundary conditions . . 102 Wave propagation in a soft PVC phantom with E = 17kP a and an axial length of 4cm under two different boundary conditions when a step excitation is applied. (a) and (c) show the displacements of tissue at different depths for a flat and a point-source excitation respectively. The location of the wavefront peaks are shown versus time for these two boundary conditions in (b) and (d), where a line has been fitted to the data to estimate the speed. The resonance frequency of the phantom is also a function of the elastic modulus, geometry and boundary conditions. . . . . . . . . . . . . . . . . . . 103 Acquisition time-line for two frames of B-mode data on a twelve element probe. . . . . . . . . . . . . . . . . . . . . . . 110 Acquisition time-line for two frames of high frame rate acquisition data on a twelve element probe . . . . . . . . . . . . . 111 A hypothetical acquisition time-line for delay-free simultaneous sampling on a twelve element probe . . . . . . . . . . . . 113 Flow chart of sequence programming for high frame rate image acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Acquisition time-line involving reprogramming delay on a twelve element probe . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 The extension of a limited duration signal to a periodic signal 120 Signal processing pipeline of the high frame rate ultrasound system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 Schematics of (a) the side view and (b) the front view of the experimental setup developed for testing the system. Note the choice of the coordinate axes. . . . . . . . . . . . . . . . . 123 Different steps in the phase correction process: the phase (a) before any compensation (b) after compensation for the delays on each scan-line (c) after compensation for the delays between the scan-lines in each sector (d) after compensation for the delays between the sectors. . . . . . . . . . . . . . . . 125 x List of Figures 5.10 (a), (b), (c), (d) the amplitude, and (e), (f), (g), (h) corresponding phase of the displacement phasors at different frequencies of excitation. From left to right the frequencies are 75, 100, 150, and 200Hz. . . . . . . . . . . . . . . . . . . . . . 5.11 Spatial derivatives of the axial displacement phasor U (x, y) at the excitation frequency of 100Hz. (a) |∂U/∂x| (b) |∂ 2 U/∂x2 | (c) |∂U/∂y| (d) |∂ 2 U/∂y 2 | (e) ∠∂U/∂x (f) ∠∂ 2 U/∂x2 (g) ∠∂U/∂y (h) ∠∂ 2 U/∂y 2 . . . . . . . . . . . . . . . . . . . . . . 5.12 Outline of the geometry of the inclusion phantoms used in the elastacity experiments. . . . . . . . . . . . . . . . . . . . . 5.13 (a) B-mode and (b) elastogram images of a phantom with a circular hard inclusion. . . . . . . . . . . . . . . . . . . . . . . 5.14 (a) B-mode and (b) elastogram images of a phantom with a circular soft inclusion. . . . . . . . . . . . . . . . . . . . . . . 6.1 126 127 128 129 129 Obtained Young’s modulus as a function of the volume portion of Plastic Softener or Hardener added during the fabrication of the phantoms . . . . . . . . . . . . . . . . . . . . . . 137 A.1 Simple schematic of the rheometry device in a conventional rheometry configuration. In this configuration the force sensor measures the force on the moving plate and moves with it. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 A.2 Simple schematic of the rheometry device in a conventional rheometry configuration. In this configuration the force sensor measures the force on the static plate and does not move. 145 A.3 The infinitesimal volume which is used in the definition of Elastic modulus. . . . . . . . . . . . . . . . . . . . . . . . . . 146 xi List of Abbreviations 1D 2D 3D A-line ARF ARFI B-mode DC motor EC ECG EDTA FEM MRE MRI MSD PC PDE PRF PVA PVA-C PVC P-wave RF ROI S-wave SNR TMM One Dimensional Two Dimensional Three Dimensional Amplitude line, the amplitude of the Hilbert transform of an RF-line Acoustic Radiation Force Acoustic Radiation Force Imaging Brightness mode, the most common type of medical ultrasound images Direct Current motor European Commission EchoCardioGram EthyleneDiamineTetraAcetic acid Finite Element Modeling or Method Magnetic Resonance Elastography Magnetic Resonance Imaging Mass Spring Damper model Personal Computer Partial Differential Equation Pulse Repetition Frequency PolyVinyl Alcohol PoylyVinyl Alcohol Cryogel PolyVinyl Chloride Primary wave, defined in seismology Radio Frequency, RF-data or RF-line is the signal received by an ultrasound transducer Region Of Interest Secondary wave, defined in seismology Signal to Noise Ratio Tissue Mimicking Material xii Acknowledgments This work would not have been possible without the consistent support of my parents, whose caring love has always been a source of motivation for me. I would also like to express my deep gratitude to my supervisors, Dr. Tim Salcudean and Dr. Robert Rohling, for giving me the freedom, yet the guidance to pursue this research to this point. I would also like to thank the members of the Robotics and Control Lab of University of British Columbia Julian Guerrero, Orcun Goksel, Ehsan Dehghan, Hani Eskandari, Reza Zahiri Azar, Sara Khosravi, and Jing Xiang for providing a collaborative and encouraging research environment. Living far away from my family, my close friends Alireza Forghani, Soraya Kasnavi, and Tricia Pang, have been the greatest source of support for me. In the end I should thank my trusted best friends Thomas Diego Prananta, Denis Tran, Andre Lei, Budi Kartono, and Joanna Lam who have supported me with their invaluable friendship throughout these years. I wish all of them the best in their lives and careers. xiii Statement of Co-Authorship This thesis is based on several manuscripts, resulting from collaboration between multiple researchers. A version of Chapter 2 appeared in the IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control in July, 2009. This article was co-authored with Dr. Hani Eskandari, Dr. Tim Salcudean, and Dr. Robert Rohling. The author’s contribution in that article was developing the main idea, numerical simulation, phantom fabrication, experimental evaluation of the algorithms and writing the manuscript. Dr. Eskandari contributed to the article through discussions, calibration of the experimental set-up, software development to control the rheometer, and editing the manuscript. Dr. Salcudean designed the rheometer platform. Dr. Rohling and Dr. Salcudean helped with their numerous suggestions in the course of the data acquisition and analysis. They also assisted by editing the manuscript. A version of Chapter 3 appeared in the Journal of Acoustical Society of America in September, 2009. This article was co-authored with Dr. Tim Salcudean and, Dr. Robert Rohling. The author’s contribution in the article was developing the idea, mathematical analysis and proof, numerical simulation, and writing the manuscript. Dr. Salcudean and Dr. Rohling assisted by editing the manuscript. A version of Chapter 4 appeared in the “Physics in Medicine and Biology” in July, 2009. This article was co-authored with Dr. Hani Eskandari, Dr. Tim Salcudean, and Dr. Robert Rohling. The original idea that the longitudinal waves can travel at a low speed was developed by the author. Dr. Eskandari made several contributions in particular adding viscosity to the analysis, analyzing the solution for a point source placed on a cylinder, and devising one of the methods for determining the wave speed from displacement spectra. The second method was the idea of the author. Dr. Eskandari also performed the FEM simulations and the analysis of the results. The experiments were first performed by the author and suggested to Dr. Eskandari. However, the actual data used in the article were collected by Dr. Eskandari. The manuscript was written by Dr. Eskandari. The author was involved at different stages of the manuscript through discussions xiv Statement of Co-Authorship with Dr. Eskandari. Dr. Salcudean and Dr. Rohling assisted with their suggestions and editing the manuscript. A version of Chapter 5 has been submitted for publication. This article was co-authored with Alex Brant, Dr. Tim Salcudean and Dr. Robert Rohling. The author’s contribution in that article was developing the main idea, and designing and carrying out the experiments. Alex Brant converted the author’s Matlab code into C++ and incorporated it into a larger program running on the ultrasound machine for real-time data acquisition. Dr. Salcudean and Dr. Rohling assisted with their suggestions throughout the course of the research and editing the manuscript. xv Chapter 1 Introduction 1.1 Motivation Each of the medical imaging modalities provides the clinicians with some type of information about the internal characteristics of tissue. The information is provided to the clinicians by the contrast in the produced images. The contrast in X-ray images is a function of the absorption and scattering properties of the tissue in the X-ray band of the electromagnetic wave spectrum. The contrast in magnetic resonance (MRI) images depends on the characteristics of the hydrogen atom in the molecules in tissue, or its concentration. The contrast in ultrasound images is a function of the absorption and scattering properties of the tissue for mechanical waves in the 1-10MHz range of frequencies. Elasticity imaging or elastography refers to medical imaging techniques which provide the clinicians with images and measurements of the viscoelastic properties of tissue [1–8]. Changes in the viscoelastic properties of tissue are often associated with pathology. In the diagnosis of breast tumors, the physician often examines the breast by palpation in order to feel the presence of hard lumps. The same is done for liver tumors or through digital rectal examination for prostate cancer. As useful as it is, palpation is not an accurate way of determining the tissue stiffness. Images of tissue elasticity, however, provide quantitative measures for the physician to make further decisions for the diagnosis or treatment. 1.2 Background In spite of their differences, all the elastography techniques developed so far share a general methodology for obtaining tissue mechanical properties. In this methodology, the tissue is subjected to some form of deformation [9–17], as it is being imaged by a medical imaging system. The images obtained from the deforming tissue are further processed to measure tissue motion, namely displacements or velocities [11, 18–27]. The tissue motion is then 1 Chapter 1. Introduction used in an inverse problem to find its mechanical properties [11,28–40]. The images of tissue elasticity are also called elastograms in the literature. When the mechanical properties of a sample are known, and the sample is subjected to forces, the motion of the material can be found by solving a forward problem. Therefore the reverse of this process, that is of finding the mechanical properties of the material from its displacements is called solving the inverse problem. The elastography techniques can be categorized based on the excitation technique, the medical imaging technique, and the inverse problem approach. Type of Excitation The excitation can be static, or dynamic. A dynamic excitation can be single frequency, or wideband. The excitation can be external, for instance a vibrator placed on the tissue surface, or induced internally with a device, or it can be natural due to heart beat or respiratory motion. [9–17]. Imaging Modality Most often, the medical imaging technique used for elastography is either ultrasound or magnetic resonance imaging (MRI). Each of these techniques has its own advantages and disadvantages. MRI can provide measurements of tissue displacement with sub-micron accuracy and a high signal to noise ratio (SNR) [11, 36, 37, 41, 42]. It can measure displacements in all three spatial directions (x, y, and z), hence provides 3D vector displacements. Moreover the 3D vector displacements can be measured at all spatial locations of a 3D volume. This is the most complete set of information about the tissue displacement, as it completely describes the way the tissue is deforming. However obtaining such MRI images can take from a few minutes to a portion of an hour. The motion of the patient during the prolonged acquisition interval reduces the SNR and the accuracy of the measurements. Moreover MRI is not as readily available as ultrasound and is more expensive to perform. Ultrasound imaging devices on the other hand are relatively inexpensive and more readily available [18–27]. The acquisition process is much faster, and with the improvements in the elastography techniques, real-time ultrasound elastography is not a far-fetched goal anymore. The disadvantages are that the measurements are less accurate than MRI and suffer from a lower SNR. The most important shortcoming of ultrasound, however, is that in its 2 Chapter 1. Introduction most basic form, it only measures 1D displacements in the direction of the ultrasound wave propagation, and the measurements are carried out only over a 2D plane, namely the plane of imaging. Even if 3D measurements could be made available through 2D or mechanical 3D arrays, ultrasound resolution in the lateral direction has fundamental limits which are typically an order of magnitude less than the axial resolution. In this case, a complete description of the tissue motion is not available and the mechanical properties need to be estimated based on the partial information. Inverse Problem Techniques The inverse problem techniques used in elastography are quite versatile and to some extent, depend on the excitation technique used. One of the simplest techniques used in manual ultrasound elastography is to show the inverse of the strain. In this technique the physician applies pressure on the tissue using the ultrasound probe itself, and it is assumed that this pressure causes a uniform constant stress in the material. Under these conditions the strain becomes inversely proportional to the tissue elasticity. Therefore in this case, the image of the inverse of the strain is in effect an elastogram. One class of the inversion techniques is based on finite element modeling (FEM) [30–34, 43–45]. In the most common approach of these techniques, the region of interest (ROI) is discretized by a mesh. An initial distribution of elasticity values is assumed for the nodes of the mesh. This could be a uniform elasticity distribution for instance. The forward problem is then solved based on these elasticity values, using FEM techniques, and the expected displacements are found. The expected displacements are then compared to the measured displacements. As these would not match, an optimization scheme is carried out to adjust the elasticity values of the nodes iteratively, in order to minimize the difference. Some of the approaches to the inverse problem are based on the mechanical wave equation [11, 13–16, 23, 36–38, 46–50]. When a dynamic excitation is used, the motion of the tissue is governed by a wave equation. The form of the equation depends on the level of modeling used for the tissue, for instance whether anisotropy or viscosity are taken into account or not. For an isotropic purely elastic homogeneous medium, in the absence of body 3 Chapter 1. Introduction forces, the equations take the form: ∂2u ∂t2 ∂2v ρ 2 ∂t ∂2w ρ 2 ∂t ρ ∂∆ + µ∇2 u ∂x ∂∆ = (λ + µ) + µ∇2 v ∂y ∂∆ = (λ + µ) + µ∇2 w ∂z = (λ + µ) (1.1) (1.2) (1.3) The wave equation is a partial differential equation; the independent variables are the spatial coordinates (x, y, z), and the time t. The dependent variables are the 3D displacement components (u, v, w). The mechanical properties of the tissue appear as coefficients: ρ is the density, and λ and µ are the Lam´e constants and describe the elastic behavior of the material. µ is also called the shear modulus of elasticity. ∆ is called the dilatation and is defined in terms of the displacements by, ∆= ∂u ∂v ∂w + + ∂x ∂y ∂z (1.4) From a control systems perspective, the inverse problem can be viewed as an identification problem; the measured quantities are the dependent variables of the differential equation, and the goal is to find the coefficients. Some of the published methods have taken this approach [51, 52]. When the 3D vector displacement field is available over a 3D volume, it is possible to carry out a full inversion of the wave equation [41, 42]. However, as is the case in ultrasound elastography, the 3D vector displacement is not readily available. Moreover, in some magnetic resonance elastography (MRE) techniques, to decrease the imaging time, not all three components of the displacements are measured. In these cases, the inversion techniques use the partial data, for instance a single component of the displacement u, measured over a plane (x, y), to estimate the mechanical properties. Most of the inversion techniques based on partial data first estimate the phase speed, then estimate the tissue elasticity [13–16, 36–38,47, 48, 50]. The phase speed is the speed at which equiphase surfaces propagate through the medium. Under certain excitation and boundary conditions, the phase speed cphase for shear waves is related to the shear modulus µ, cphase = µ . ρ (1.5) 4 Chapter 1. Introduction Different techniques exist for estimating the phase speed. In transient elastography, the excitation is in the form of a short pulse [13–15]. The propagation speed of this pulse is then measured. In some of these techniques the measured quantity is exactly the phase speed. In another class of dynamic elastography techniques, the excitation is continuously applied to the tissue, and the measured displacements are in the steady state [11, 36, 37, 48]. In these techniques the displacement are analyzed in the phasor domain. Various techniques such as algebraic inversion and phase gradient have been proposed for estimating the phase speed. The underlying problem which all these techniques strive to solve is that of estimating the local frequency of a sampled signal with a finite duration. This problem is an inherently challenging problem [53, 54]. As mentioned before, when the 3D vector displacements are available over a 3D volume, a complete inversion can be carried out. However in the case of partial displacement data, where the phase speed is estimated, a relationship exists between the resolution of the elastograms and the frequency of excitation. This is due to the fact that to measure the local frequency of a signal with a higher accuracy, more of its periods should be observed. As a result, the higher the frequency of excitation, the shorter the wavelength of the waves, and the higher the resolution of the elastograms. Elastograms with higher resolution are clinically preferred, because they can show smaller tumors and calcifications and lead to earlier diagnosis. But the higher frequencies are attenuated more in real tissue, so there is a practical limit to the number of wavelengths that can be observed in a given ROI. Nevertheless, there are still limitations on the ability of the ultrasound machine to capture motion at these frequencies. 1.3 Thesis Objectives This thesis studies the inversion problem of elastography under the following conditions: • Excitations are continually applied to the tissue. • The displacements can be analyzed in the phasor domain. • The imaging system is an ultrasound system. • Only partial measurements of displacements are carried out. • The phase speed is used for estimating the tissue elasticity. The objective of this thesis is to establish the hypothesis that a wave equation-based approach can be used in ultrasound elastography to provide images of tissue elasticity. To examine this hypothesis, the following 5 Chapter 1. Introduction objectives are set: • Finding a method for accurately measuring the elasticity of tissue and tissue mimicking material samples for setting gold standards. • Fabricating tissue mimicking materials with controllable elasticity. • Establishing the relationship between the phase speed and tissue elasticity through a theoretical analysis of the elastic wave equation. • Developing an ultrasound elastography system which uses the developed algorithm for forming elastograms. The following contributions have been made towards achieving the objectives: • A novel rheometry method has been devised which is capable of accurate measurements of tissue elasticity and viscosity on test blocks. The accuracy of the device has been shown by a comparison of the results with conventional rheometry. • Poly Vinyl Chloride (PVC) tissue mimicking materials have been fabricated and their elasticity characterized based on the volume portion of added softener or hardener. • The limitations of the shear wave speed equation (1.5) have been established. In particular experimental data have been provided for which equation (1.5) does not apply. • It has been established theoretically that a universal relationship does not exist between the phase speed and tissue elasticity. However under special boundary and excitations, a relationship can exist. • The relationship between the phase speed and tissue elasticity is established for small and large profile exciters placed on finite media under different boundary conditions, by using FEM simulations and also experimentally. • A novel high-frame-rate ultrasound system has been devised which uses the idea of sector subdivision to achieve the increased frame rate. Unlike other systems that use sector subdivision, the devised system compensates for the time delays introduced in the acquisition process by a novel frequency domain phase compensation scheme. • Conditions have been found which should be met by the mechanical excitation in order for the high-frame-rate system to work properly. • The high-frame-rate system has been used in conjunction with a phase speed based inversion algorithm to develop and ultrasound elastography system. The performance of the elastography system has been tested on tissue mimicking material. 6 Chapter 1. Introduction 1.4 Chapter Summary The overall goal of this thesis is to study the application of wave equationbased inversion techniques in ultrasound elastography. This thesis has been written in the manuscript-based style, as permitted by the Faculty of Graduate Studies at the University of British Columbia. In the manuscript-based thesis, each chapter represents an individual work that has been published, submitted or prepared for submission to a peer reviewed publication. Each chapter is self-sustained in the sense that it includes an introduction to the work presented in that chapter, the methodology, results and discussion. The references are summarized in the bibliography found at the end of each chapter. The appendices pertaining to each chapter are presented at the end of the thesis. In Chapter 2 we show that for a finite block of material, the phase speed for longitudinal waves could be related to the shear modulus as in, cphase = µ 3 . ρ (1.6) Furthermore, based on this relationship, a novel method is introduced for rheometry, i.e. determining the elasticity and viscosity of sample blocks of material. This method is more suitable for studying soft tissue specimens than conventional rheometry. The accuracy of the method is tested against conventional rheometry. In Chapter 3 a further theoretical analysis is carried out on the phase speed which reveals that neither (1.5) nor (1.6) is a universal relation. In fact, for the same material, with the same shear modulus, the phase speed can take different values, depending on the beam form of the wave, and the boundary conditions. In Chapter 4 the results of Chapter 3 are corroborated by using FEM and experiments. With a 3D FEM simulation, it is shown that a point exciter creates longitudinal waves which travel at the phase speed given by (1.5) but a large profile exciter creates longitudinal waves which travel at the phase speed given by (1.6). Both harmonic and transient analyses are performed. Similar results are obtained from phantom experiments. In Chapter 5 a method is presented for acquiring displacement phasors at a high-frame-rate. The low frame rate of conventional ultrasound systems limits the scope of motions they can track. Namely, the Nyquist criterion imposes a limit on the highest frequency of excitation that can be used with these systems. This limit is in the range of 50Hz to 200Hz, depending on 7 Chapter 1. Introduction the depth of imaging and field of view. In the current setting, to form elastograms of higher resolution it is necessary to excite the tissue at higher frequencies, which in turn requires a higher frame rate on the imaging system. The proposed method can easily be implemented on most commercial ultrasound systems without the need for a hardware upgrade. The effectiveness of the system in forming elastograms is demonstrated experimentally on tissue mimicking material with soft and hard inclusions. In Chapter 6 the thesis is concluded by a brief summary and a discussion which shows the place of the contributions of this thesis in the bigger picture of tissue elastography. Recommendations are made for future research directions. Two appendices are appended to the thesis. In Appendix A the limitations of conventional rheometry are analyzed. In Appendix B a literature survey of the most commonly used tissue mimicking materials is presented together with a section on fabricating Poly Vinyl Chloride phantoms. In Appendix C the method for deriving the phasor of the displacements is described together with two typical inversion techniques based on the wave equation. 8 References [1] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: A quantitative method for imaging of elasticity of biological tissues,” Ultrasonic Imaging, vol. 13, pp. 111–134, 1991. [2] K. Parker, L. Gao, R. Lerner, and S. Levinson, “Techniques for elastic imaging: A review,” IEEE Engineering in Medicine and Biology Magazine, vol. 15, no. 6, pp. 52–59, 1996. [3] L. Gao, K. J. Parker, R. M. Lerner, and S. Levinson, “Imaging of the elastic properties of tissue-a review,” Ultrasound Med. Biol., vol. 22, no. 8, pp. 957–977, 1996. [4] J. Ophir, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, R. Righetti, and T. Varghese, “Elastographic imaging,” Ultrasound Med. Biol., vol. 26, pp. S23–S29, 2000. [5] A. Sarvazyan, Handbook of Elastic Properties of Solids, Liquids and Gases. Academic Press, 2001, vol. III, ch. Elastic Properties of Soft Tissues. [6] J. Ophir, S. K. Alam, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, C. Merritt, R. Righetti, R. Souchon, S. Srinivasan, and T. Varghese, “Elastography: Imaging the elastic properties of soft tissues with ultrasound,” Journal of Medical Ultrasonics, vol. 29, pp. 155–171, 2002. [7] J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imaging elastic properties of biological tissues,” Annual Review Biomedical Engineering, vol. 5, pp. 57–78, April 2003. [8] K. J. Parker, L. S. Taylor, S. Gracewski, and D. J. Rubens, “A unified view of imaging the elastic properties of tissue,” J. Acoust. Soc. Am., vol. 117, no. 5, pp. 2705–2712, May 2005. [9] R. Lerner and K. Parker, “Sonoelasticity images, ultrasonic tissue characterization and echographic imaging,” in Proceedings of the 7th European Communities Workshop, 1987, nijmegen, Netherlands. 9 Chapter 1. Introduction [10] R. Lerner, K. Parker, J. Holen, R. Gramiak, and R. Waag, “Sonoelasticity: Medical elasticity images derived from ultrasound signals in mechanically vibrated targets,” Acoustical Imaging, vol. 16, pp. 317– 327, 1988. [11] R. Sinkus, J. Lorenzen, D. S. amd M. Lorenzen, M. Dargatz, and D. Holz, “High-resolution tensor mr elastography for breast tumor detection,” Physics Med. Biol., vol. 45, pp. 1649–1664, 2000. [12] R. Souchon, L. Soualmi, M. Bertrand, J.-Y. Chapelon, F. Kallel, and J. Ophir, “Ultrasonic elastography using sector scan imaging and a radial compression,” Ultrasonics, vol. 40, pp. 867–871, 2002. [13] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulus imaging with 2-d transient elastography,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 49, no. 4, pp. 426–435, April 2002. [14] L. Sandrin, M. Tanter, J.-L. Gennisson, S. Catheline, and M. Fink, “Shear elasticity probe for soft tissues with 1-d transient elastography,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 49, no. 4, pp. 436– 446, April 2002. [15] J.-L. Gennisson, S. Catheline, S. Chaffa¨ı, and M. Fink, “Transient elastography in anisotropic medium: Application to the measurement of slow and fast shear wave speeds in muscles,” J. Acoust. Soc. Am., vol. 114, no. 1, pp. 536–541, July 2003. [16] J. Bercoff, R. Sinkus, M. Tanter, and M. Fink, “Supersonic shear imaging: A new technique for soft tissue elasticity mapping,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 51, no. 4, pp. 396–409, 2004. [17] M. W. Urban and J. F. Greenleaf, “Harmonic pulsed excitation and motion detection of a vibrating reflective target,” J. Acoust. Soc. Am., vol. 123, no. 1, pp. 519–533, Jan 2008. [18] F. Yeung, S. F. Levinson, D. Fu, and K. J. Parker, “Feature-adaptive motion tracking of ultrasound image sequences using a deformable mesh,” IEEE Trans. Med. Imaging, vol. 17, no. 6, pp. 945–956, 1998. [19] E. Konofagou, T. Varghese, J. Ophir, and S. Alam, “Power spectral strain estimators in elastography,” Ultrasound Med. Biol., vol. 25, no. 7, pp. 1115–1129, 1999. 10 Chapter 1. Introduction [20] E. E. Konofagou and J. Ophir, “Precision estimation and imaging of normal and shear components of the 3d strain tensor in elastography,” Phys. Med. Biol., vol. 45, pp. 1553–1563, 2000. [21] T. Varghese, E. Konofagou, J. Ophir, S. Alam, and M. Bilgen, “Direct strain estimation in elastography using spectral cross-correlation,” Ultrasound Med. Biol., vol. 26, no. 9, pp. 1525–1537, 2000. [22] M. Fink, L. Sandrin, M. Tanter, S. Catheline, S. Chaffai, J. Bercoff, and J. Gennisson, “Ultra high speed imaging of elasticity,” in IEEE Ultrasonics Symposium, 2002, pp. 1811–1820. [23] J. Bercoff, R. Sinkus, M. Tanter, and M. FinkJ, “3d ultrasound-based dynamic and transient elastography: First in vitro results,” in IEEE Ultrasonics Symposium, 2004, pp. 28–31. [24] K. Hoyt, F. Forsberg, and J. Ophir, “Investigation of parametric spectral estimation techniques for elasticity imaging,” Ultrasound Med. Biol., vol. 31, no. 8, pp. 1109–1121, 2005. [25] ——, “Comparison of shift estimation strategies in spectral elastography,” Ultrasonics, vol. 44, pp. 99–108, 2006. [26] ——, “Analysis of a hybrid spectral strain estimation technique in elastography,” Phys. Med. Biol. 51, vol. 51, pp. 197–209, 2006. [27] J.-L. Gennisson, T. Deffieux, R. Sinkus, P. Annic, M. Pernot, F. Cudeiro, G. Montaldo, M. Tanter, M. Fink, and J. Bercoff, “A 3d elastography system based on the concept of ultrasound-computed tomography for in vivo breast examination,” in IEEE Ultrasonics Symposium, 2006, pp. 1037–1040. [28] A. J. Romano, J. J. Shirron, and J. A. Bucaro, “On the noninvasive determination of material parameters from a knowledge of elastic displacements: Theory and numerical simulation,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 45, no. 3, pp. 751–759, May 1998. [29] C. Sumi, A. Suzuki, and K. Nakayama, “Estimation of shear modulus distribution in soft tissue from strain distribution,” IEEE Trans. Biomed. Eng., vol. 42, no. 2, pp. 193–202, February 1995. [30] K. Raghavan and A. E. Yagle, “Forward and inverse problems in elasticity imaging of soft tissues,” IEEE Transactions on Nuclear Science, vol. 41, no. 4, pp. 1639–1648, August 1994. 11 Chapter 1. Introduction [31] F. Kallel and M. Bertrand, “Tissue elasticity reconstruction using linear perturbation method,” IEEE Trans. Med. Imaging, vol. 15, no. 3, pp. 299–313, June 1996. [32] M. Doyley, P. Meaney, and J. Bamber, “Evaluation of an iterative reconstruction method for quantitative elastography,” Phys. Med. Biol., vol. 45, pp. 1521–1540, 2000. [33] P. E. Barbone and J. C. Bamber, “Quantitative elasticity imaging: What can and cannot be inferred from strain images,” Phys. Med. Biol., vol. 47, pp. 2147–2164, 2002. [34] A. A. Oberai, N. H. Gokhale, and G. R. Feij´oo, “Solution of inverse problems in elasticity imaging using the adjoint method,” Inverse Problems, vol. 19, pp. 297–313, 2003. [35] D. Fu, S. Levinson, S. Gracewski, and K. Parker, “Non-invasive quantitative reconstruction of tissue elasticity using an iterative forward approach,” Phys. Med. Biol. 45 (2000) 14951509, vol. 45, pp. 1495–1509, 2000. [36] A. Manduca, T. Oliphant, M. Dresner, J. Mahowald, S. Kruse, E. Amromin, J. Felmlee, J. Greenleaf, and R. Ehman, “Magnetic resonance elastography: Non-invasive mapping of tissue elasticity,” Medical Image Analysis, vol. 5, pp. 237–2540, 2001. [37] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complexvalued stiffness reconstruction by for magnetic resonance elastography by algebraic inversion of the differential equation,” Magnetic Resonance in Medicine, vol. 45, pp. 299–310, 2001. [38] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelkaram, and J. Culiolic, “Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: An inverse problem approach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec 2004. [39] B. Robert, R. Sinkus, B. Larrat, M. Tanter, and M. Fink, “A new rheological model based on fractional derivatives for biological tissues,” in IEEE Ultrasonics Symposium, 2006, pp. 1033–1036. [40] X. Zhang and J. F. Greenleaf, “Estimation of tissues elasticity with surface wave speed (l),” J. Acoust. Soc. Am., vol. 122, no. 5, pp. 2522– 2525, November 2007. 12 Chapter 1. Introduction [41] R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink, “Viscoelastic shear properties of in vivo breast lesions measured by mr elastography,” Magnetic Resonance Imaging, vol. 23, pp. 159–165, 2005. [42] R. Sinkus, M. Tanter, S. Catheline, J. Lorenzen, C. Kuhl, E. Sondermann, and M. Fink, “Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography,” Magnetic Resonance in Medicine, vol. 53, no. 2, pp. 372–387, 2005. [43] A. A. Oberai, N. H. Gokhale, M. M. Doyley, and J. C. Bamber, “Evaluation of the adjoint equation based algorithm for elasticity imaging,” Phys. Med. Biol., vol. 49, pp. 2955–2974, 2004. [44] M. M. Doyley, S. Srinivasan, S. A. Pendergrass, Z. Wu, and J. Ophir, “Comparative evaluation of strain-based and model-based modulus elastography,” Ultrasound Med. Biol., vol. 31, no. 6, pp. 787–802, 2005. [45] M. M. Doyley, S. Srinivasan, E. Dimidenko, N. Soni, and J. Ophir, “Enhancing the performance of model-based elastography by incorporating additional a priori information in the modulus image reconstruction process,” Phys. Med. Biol., vol. 51, pp. 95–112, 2006. [46] J. Bercoff, S. Chaffai, M. Tanter, L. Sandrin, S. Catheline, M. Fink, J. L. Gennisson, and M. Meunier, “In vivo breast tumor detection using transient elastography,” Magnetic Resonance in Medicine, vol. 29, no. 10, pp. 1387–1396, 2003. [47] J. Bercoff, M. Muller, M. Tanter, and M. Fink, “Study of viscous and elastic properties of soft tissue using supersonic shear imaging,” in IEEE Ultrasonics Symposium, 2003, pp. 925–928. [48] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity and viscosity from measurement of shear wave speed dispersion,” J. Acoust. Soc. Am., vol. 115, no. 6, pp. 2781–2785, June 2004. [49] J. Fromageau, J.-L. Gennisson, C. Schmitt, R. L. Maurice, R. Mongrain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel mechanical properties with four ultrasound elastography methods and comparison with gold standard testing,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007. [50] A. Romano, J. Bucaro, P. Abraham, and S. Dey, “Inversion methods for the detection and localization of inclusions in structures utilizing 13 Chapter 1. Introduction dynamic surface displacements,” in Proceedings of SPIE Vol. 5503, vol. 5503, 2004, pp. 367–374. [51] E. Turgay, S. Salcudean, and R. Rohling, “Identifying mechanical properties of tissue by ultrasound,” Ultrasound Med. Biol., vol. 32, no. 2, pp. 221–235, 2006. [52] H. Eskandari, S. Salcudean, and R. Rohling, “Viscoelastic parameter estimation based on spectral analysis,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 55, no. 7, pp. 1611–1625, 2008. [53] B. Boashash, “Estimating and interpreting the instantaneous frequency of a signal – part 1:fundamentals,” Proceedings of IEEE, vol. 80, no. 4, pp. 520–538, April 1992. [54] ——, “Estimating and interpreting the instantaneous frequency of a signal – part 1:fundamentals,” Proceedings of IEEE, vol. 80, no. 4, pp. 540–568, April 1992. 14 Chapter 2 Measurement of Viscoelastic Properties of Tissue Mimicking Material Using Longitudinal Wave Excitation∗ 2.1 Introduction The viscoelastic characterization of tissue and tissue mimicking material is of interest in many areas of biomedical research. The variability in tissue stiffness due to anatomical features, as well as the changes incurred in tissue stiffness due to pathology, has led to the development of elastography as a clinical imaging modality [1–5]. In general, elastography involves observation of mechanical deformation patterns inside tissue by means of a clinical imaging system and use of these patterns to infer mechanical properties of the tissue being imaged. To create deformation inside tissue different methods have been devised. Acoustic excitation by sound waves was used by Lerner et al. [6,7]. Ophir et al. proposed quasi-static compression [8]. Rudenko et al. used acoustic radiation force of focused ultrasound to create shear waves inside tissue [9,10]. Sandrin et al. used external shakers which create longitudinally polarized shear waves [11–13]. Excitation by two or more sources of slightly different frequencies has also been considered. The sources are focused ultrasound ∗ A version of this chapter has been published: A. Baghani, H. Eskandari, S.E. Salcudean, R. Rohling, “Measurement of Viscoelastic Properties of Tissue Mimicking Material Using Longitudinal Wave Excitation,” IEEE Trans. on Ultras. Ferr. Freq. Cont., vol. 56, no. 7, pp. 1405-1418, Jul. 2009 15 Chapter 2. Elastography Using Longitudinal Waves in the work by Fatemi and Greenleaf [14]. They are shear wave sources in the work by Hoyt et al. [15]. Among these works, only Ophir et al. used compressional excitation, and only in the quasi-static settings. In all the other cases, either dominantly shear excitation is used, or the compressional excitation is ignored or suppressed by taking the curl of the displacement field. To quantitatively study the effectiveness of elastography imaging methods in providing images of viscoelastic properties of tissue, it is of interest to test these methods on phantoms with known anatomical features of known viscoelastic properties. A number of researchers are involved in the development of phantoms with adjustable viscoelastic properties [16–23]. Some researchers have used conventional rheometry to validate the results obtained from their elastography techniques [24,25]. Other researchers have developed specific methods for measuring viscoelastic properties of the phantom and soft tissue specimens [26, 27]. Srinivasan et al. used a nanoindenter to measure elasticity [28] . Kalanovic et al. compared the results of elasticity measurements with a commercial indenter to the results obtained from their developed shear rheometry device [29]. Egorov et al. have developed a soft tissue elastometer which uses a small profile indenter and a scale to measure the applied force [30]. They concluded that the validity of their results depends on the size of the tissue sample relative to the size of the indenter. Samani et al. also used a small profile indenter to measure elasticity [31]. FEM was used in their work to study the geometric constant relating the force to the indentation. In their work magnetic resonance imaging based elastography was also proposed for the purpose of rheometry. Mechanical wave propagation in elastic continua has been studied extensively and is well understood [32, 33]. Some viscoelastic continua have also been studied with certain assumptions about the viscoelastic behavior of the matter [34,35]. In general the stress-strain relationship for a material can be very complex. The assumption usually made is that the stress-strain relationship can be modeled by linear time invariant differential equations. Low frequency shear waves were used in the work of Catheline et al. for the measurement of the viscoelastic properties of tissue [36]. Longitudinally polarized transient shear waves were produced mechanically and the resulting displacements were tracked by ultrasound speckle tracking to give measurements of the wave speed and attenuation through an inverse problem approach. It was shown that a Voigt model fits the visoelastic behavior of the living tissue better than a Maxwell model. It is the objective of this research to: 16 Chapter 2. Elastography Using Longitudinal Waves • Show theoretically and experimentally that low frequency longitudinal waves can propagate in a tissue mimicking material of finite dimensions at speeds which are much lower than those normally associated with ultrasound propagation. • Show experimentally that it is possible to use ultrasound-based speckle tracking techniques to track the longitudinal waves. • Show theoretically and experimentally that the measured displacements caused by longitudinal wave propagation can be used to estimate the viscoelastic properties of the tissue mimicking material. To this end, a new framework is developed for viscoelastic characterization of tissue and phantom specimens by using longitudinal wave excitation and validated against the conventional rheometry. We start by a brief review of longitudinal wave equations in elastic media in Section 2.2. The experimental setup, material, and methods used are presented in Section 2.3. Section 2.4 is devoted to the presentation of our rheometry methods. The experimental results are presented in Section 2.5 and discussed in Section 2.6. We conclude the paper in Section 2.7. 2.2 Theory of Wave Motion In this section the theory of wave motion in elastic media is briefly reviewed and the notation used in the article presented. The presentation is adopted from [33]. We start with the model of an elastic medium. We use u(x) to denote the displacement vector field. τij and ij are used to denote the stress and strain tensors. By definition: ij = 1 2 ∂ ∂ ui + uj ∂xj ∂xi . (2.1) The fundamental conservation equations in a medium can then be expressed as: mass : momentum : angular momentum : energy : ∂ ∂t ρ + ∂ ˙ i) i ∂xi (ρu ∂ j ∂xj τji = 0, (2.2) + ρfi = ρ¨ ui , (2.3) τij = τji ρe˙ = i=j, ij τji ˙ij . (2.4) (2.5) where ρ(x) is the density, f (x) is the force per unit mass vector field, and e(x) is the internal energy per unit mass. These equations hold regardless of 17 Chapter 2. Elastography Using Longitudinal Waves the type of the medium. To model the specific behavior of different media, they are concatenated by a set of constitutive equations. These equations model the stress-strain behavior of the medium. If we assume a homogeneous isotropic elastic medium, under small strains, the constitutive equations widely used are τij = λ( 11 + 22 + 33 )δij + 2µ ij , (2.6) where δij is equal to one if i = j, and zero otherwise. λ and µ are called the Lam´e constants. A number of other constants are also defined based on the Lam´e constants, E = µ(3λ + 2µ)/(λ + µ) = Young s modulus , 1 ν = λ/(λ + µ) = Poisson s ratio , 2 2 K = λ + µ = bulk modulus . 3 (2.7) (2.8) (2.9) The inverse of the equation (2.6) can be expressed in terms of these constants as 1+ν ν τij − (τ11 + τ22 + τ33 )δij . (2.10) ij = E E At this point, it is worth defining the scope of our problem in terms of these parameters in dealing with soft human tissue and tissue mimicking material. It has been found experimentally that these materials are nearly incompressible [37,38]. Hence the Poisson’s ratio ν for these materials is very close to 0.5. In this respect, the material under study is closer to liquids such as water and hydraulic oils than to solids such as steel and aluminum. In terms of having a definite shape, a tissue mimicking material is closer to solids than liquids. A number of consequences follow from these peculiar properties: ν ≈ λ 0.5 , (2.11) µ, (2.12) E ≈ 3µ , (2.13) K ≈ λ. (2.14) In the rest of this article it is assumed that the relations (2.11) to (2.14) hold. 18 Chapter 2. Elastography Using Longitudinal Waves The equations (2.2) to (2.6) form a complete set of equations for the displacement fields, stresses and density and can be used to derive the most general form of the wave equation [32, 33]. For our purposes it suffices to consider some special cases. We are particularly interested in longitudinal (compressional) waves. The two cases of interest are plane waves in an infinite medium, and plane waves in a finite block of medium. 2.2.1 Longitudinal Plane Waves in an Infinite Medium In this case, it is assumed that particle motion is only in one direction (x1 direction), i.e. u2 = u3 ≡ 0. Also because of the symmetry of the plane wave, there is no variation of u1 in the x2 and x3 directions, i.e. ∂u1 /∂x2 = ∂u1 /∂x3 ≡ 0. These assumptions can be used to derive the wave equation for a longitudinal plane wave propagating in the x1 direction in an infinite elastic medium, ∂2 ∂2 (λ + 2µ) 2 u1 = ρ 2 u1 . (2.15) ∂t ∂x1 This equation predicts that the wave speed would be cinf = λ + 2µ . ρ (2.16) Because of the assumptions, in this case the strain tensor takes the form 11 ij ∂ u1 , ∂x1 = 0 otherwise . = (2.17) (2.18) The stress tensor can be found using the constitutive equations (2.6), τ11 = (λ + 2µ) 11 , (2.19) τ22 = λ 11 , (2.20) τ33 = λ 11 , (2.21) τij = 0 otherwise . (2.22) It is evident from these equations that an infinitesimal cubic volume in the medium experiences normal forces acting on all of its surfaces as the plane wave propagates through that small volume. Before we move on to the next case, it is important to mention that 19 Chapter 2. Elastography Using Longitudinal Waves the wave speed which is widely reported and accepted in the literature for ultrasound is λ + 23 µ K cUS = = (2.23) ρ ρ where K is the adiabatic compressibility in the case of fluids and bulk modulus in the case of solids [39]. If the tissue is modeled as a fluid and the longitudinal wave propagation as an adiabatic phenomenon, the wave speed would be cUS . If it is modeled as an elastic solid, the wave speed would be cinf of equation (2.16). However the difference between these two numbers is much smaller than the accuracy of the practical measurements in this research, so cinf will be used throughout this paper as the ultrasound speed. Based on ultrasound theory and a large body of evidence, when a small source of longitudinal vibrations (piezoelectric crystal) with sub-millimeter dimensions is placed on the skin, with the human body as the propagation medium, the pulses travel with a speed given by (2.23) or (2.16). Therefore for a sub-millimeter wave source of high frequency, the body acts as an infinite medium. As will be shown in the next subsection, the same is not true for all cases of longitudinal wave propagation in solids. 2.2.2 Longitudinal Plane Waves in a Finite Medium Another case of special interest is the propagation of longitudinal waves in material blocks. In this case the effects of finiteness of media should also be taken into account. More specifically, if a rectangular block of material is subjected to a spatially uniform, possibly time varying force on one of its surfaces, a longitudinal wave would propagate in the block. The study of waves of this type could be facilitated by making some assumptions about the stress and strain patterns inside the block. Two main assumptions are made here: that the “strain is uniform on each cross section of the block at every instant”, and that “plane sections remain plane under stress”. These assumptions make it possible to model the phenomenon as a one-dimensional problem. We further elaborate on these assumptions. A medium with a finite cross section acts as a waveguide for elastic waves. At each frequency of excitation a number of spatial modes is excited. Each of these modes has its own wave speed associated with it. As a result the wave speed depends on the frequency of excitation. Now consider an infinitely long rod of finite cross section which experiences a spatially uniform low frequency compression over its cross section. Under these conditions the dominant mode of the rod 20 Chapter 2. Elastography Using Longitudinal Waves has a wave pattern which satisfies the assumptions made here [32]. More precisely, at low frequencies only the first mode is excited and the amplitude of vibrations is almost constant over the cross section. This justifies the use of the assumptions made here and the use of a 1D model. Now assume that the wave is propagating in the x1 direction. The fact that the block has free boundaries on its sides causes the lateral stresses to vanish on the surfaces, τ22 = 0 , (2.24) τ33 = 0 , (2.25) τij = 0 i=j. (2.26) Because of the assumptions made, these conditions hold throughout the block. If not, the infinitesimal volumes located around the boundary of the block, would be (axially) compressed more, because of the absence of lateral stresses, than internally located infinitesimal volumes. As a result plane sections would not remain plane anymore which contradicts the assumptions. Moreover, the top and bottom boundaries of the block should be frictionless slip boundaries to prevent any external lateral forces on them. Now the constitutive equation (2.10) gives 11 = = 1+ν ν τ11 − τ11 E E 1 τ11 , E (2.27) or equivalently τ11 = E 11 µ(3λ + 2µ) = λ+µ 11 . (2.28) As a result, the equation of a longitudinal wave propagating in the x1 direction in the block would be µ(3λ + 2µ) ∂ 2 ∂2 u1 = ρ 2 u1 . 2 λ + µ ∂x1 ∂t (2.29) This equation predicts that longitudinal waves propagate in the block at a 21 Chapter 2. Elastography Using Longitudinal Waves speed of cfin = µ(3λ+2µ) λ+µ ρ = E . ρ (2.30) A comparison of (2.30) with (2.16) shows that the propagation speed of longitudinal (compressional) waves is different in finite blocks as compared to infinite media of √ the same material. For nearly incompressible material cfin ≈ 3µ/ρ = 3cS where cS = µ/ρ is the shear wave speed in the infinite medium. This speed is in the order of 1 to 10m/s for soft tissue which is much smaller than cinf which is approximately 1540m/s. Since cfin is much smaller than cinf , it is possible to track the longitudinal wave propagation in blocks of tissue using pulsed ultrasound. Our objective is to perform rheometry measurements by studying the longitudinal wave patterns inside finite blocks of tissue phantoms. In the next section the apparatus and methods used to create these wave patterns and collect data are presented. The proposed rheometry methods for analyzing the collected data follow in Section 2.4. 2.3 Experimental Methods and Setup In this section the experimental setup is discussed, and the acquired data are described. For ultrasound elastography studies, a device has been developed which is capable of applying dynamic compressions to a phantom while it is being imaged by an ultrasound probe. The detailed description of this device can be found in [40]. Figure 2.1 shows a simplified schematic of this device. The device has two plates between which the phantom is compressed. The upper plate is connected to the chassis and does not move. This plate has an opening for an ultrasound probe, and fitted by a linear array ultrasound probe. The lower plate is confined to move in one direction by a linear bearing. This plate is connected to a cam follower mechanism sitting on a cam driven by a controlled DC motor. For testing, a rectangular phantom block is placed between the two plates and the height of the upper plate is adjusted to have some pre-compression on the phantom to ensure contact with the ultrasound transducer is maintained. By (2.24) - (2.26) it is required that the top and bottom surfaces of the block only experience forces in the vertical direction. To match the experimental conditions to the theoretical ones, ultrasound gel was used as a lubricant on the top and bottom surfaces of the phantom in contact with 22 Chapter 2. Elastography Using Longitudinal Waves Figure 2.1: Simple schematic of the rheometry device the plates to reduce the friction. Super Soft PlasticTM (M-F Manufactureing Co., Inc. Fort Worth, Tx, USA), a PVC (Poly Vinyl Chloride) based plastic, was used as the phantom material. The stiffness of this material can be adjusted by adding Plastic SoftenerTM or Plastic HardenerTM provided by the same company. SigmaCell r Cellulose, Type 50 (Sigma-Aldrich Inc., St Louis, MO, USA) was added to the phantoms as the ultrasound scatterer to develop the speckle pattern needed for ultrasound-based motion tracking. For each 100mL of liquid plastic, 1mL of cellulose was added. We used a 40mm×30mm×30mm mold to fabricate five homogeneous blocks with different stiffnesses. The fabrication parameters, together with the measured densities and ultrasound speeds cinf for these phantoms are summarized in Table 2.1. The densities were measured by weighing the phantoms with a precision scale, and measuring the block dimensions with calipers and dividing the mass by the calculated volume. The ultrasound speed was found by first measuring the depth to the bottom of the phantom on a B-mode ultrasound image which is displayed on the ultrasound machine assuming an ultrasound speed of 1540m/s, and then correcting the speed so that the measured depth matches the height of the phantom. As can be seen, the density and cinf are not greatly affected by changing the proportion of Plastic to Softener or Hardener. In all the experiments the lower plate was excited with Gaussian white23 Chapter 2. Elastography Using Longitudinal Waves Table 2.1: Fabrication and measured properties of the five Plastic Hardener Softener density phantom (volume (volume (volume (kg/m3 ) ±5% No portion) portion) portion) 1 4 0 2 1020 2 5 0 1 1030 3 6 0 0 990 4 5 1 0 1010 5 4 2 0 1030 test phantoms ultrasound speed (m/s)±1% 1420 1420 1410 1440 1410 Table 2.2: Parameters of the ultrasound imaging system parameter value probe type L9–4 linear array center frequency 4MHz sampling frequency 20MHz depth of focus 25mm imaging depth 40mm Doppler pulse repetition frequency 2.5kHz number of samples acquired 75,472 noise, band-limited to 4–150Hz, while RF-lines were collected for 30 seconds at a rate of 2500 lines/second in the Doppler mode from a single A-line. The imaging system used was a Sonix RP ultrasound machine (Ultrasonix Medical Corp., Richmond, BC, Canada). The parameters of the imaging system are summarized in Table 2.2. The acquired RF-line data were processed off-line to obtain the displacement profile of each point along the scan line. Since the models used are 1D models, only the axial displacements were calculated. Lateral displacements also occur which could be used to improve the estimations in future work. The algorithm used for displacement calculations is described in [41]. The parameters used are summarized in Table 2.3. We use the notation um (x, t) to denote the measured displacement of a point located at depth x at time t. It is worth mentioning that both x and t are discrete variables in this notation. The discrete Fourier transform was then applied to um (x, t) to obtain the phasor displacement profile Um (x, jω). Frequency domain averaging was also used, as provided by the tfestimate function of MATLAB r to 24 Chapter 2. Elastography Using Longitudinal Waves Table 2.3: Parameters of the speckle tracking motion estimation algorithm parameter value number of blocks 80 block length 1mm block overlap 50% displacement measurement relative (corresponding to block velocity) stretching not used smooth the data. The data were normalized at each frequency separately so that the element located at x = 40mm would have unit amplitude and zero phase. This normalization does not have any effect on the results, but it helps with the plotting of the phasor data. At each frequency ω, Um (x, jω) is a complex-valued function of x which gives the steady-state displacement phasor of a point located at x inside the phantom, if the phantom were excited with frequency ω. A typical phasor profile is plotted in Fig. 2.2. Parts (a) and (b) show the amplitude and phase of the phasor profile as a function of depth x, while part (c) shows the same profile as a curve in the complex plane (parameterized by the depth x). In the next section Um (x, jω) will be used to characterize the viscoelastic properties of the phantoms. To validate our estimation results, we collected conventional rheometry data from the phantoms. This was achieved by a modification to the device of Fig. 2.1, where the ultrasound probe and upper plate were replaced by a force sensor (Load Cell MDB-2.5, Transducer Techniques, Temecula, CA, USA). The rheometry device was reported in a previous article [40]. More details about the rheometry method used in this article can be found in the Appendix A. 2.4 Proposed Rheometry Methods In this section two methods are introduced for analyzing the obtained data presented in the previous section for estimation of viscoelastic properties of tissue phantoms. We start by deriving a viscoelastic model for the phantom block to be used in these methods. Longitudinal wave propagation in a viscoelastic block can be modeled in a number of different ways [32, 33]. The most complete model is obtained by solving the three dimensional wave equation in a viscoelastic medium with boundary conditions adopted to match that of the block. Except for 25 Chapter 2. Elastography Using Longitudinal Waves (a) (b) (c) Figure 2.2: A typical normalized displacement phasor profile, belonging to phantom No.3 at a frequency of 90Hz. The (a) magnitude and (b) phase of the phasor profile (c) the same normalized phasor profile in the complex plane. The profile is normalized to the phasor of the point at x = 40mm. 26 Chapter 2. Elastography Using Longitudinal Waves some special cases for which analytical solutions have been derived [33], this model is quite difficult to solve and requires numerical methods such as finite element methods. Thus approximate models are usually used. The simplest approximate model is the model presented in equation (2.29) of subsection 2.2.2. From now on for brevity we denote the first coordinate by x ≡ x1 . In the frequency domain, this equation becomes (jω)2 ρU µ(3λ + 2µ) ∂ 2 U λ + µ ∂x2 ∂2 = E 2 U. ∂x = (2.31) A next step in the modeling is the introduction of a first order viscosity term. The Voigt model was shown to more closely follow the behavior of the soft tissue [36]. In this case the wave equation becomes (jω)2 ρU ∂2 U ∂x2 ∂2 +Bjω 2 U , ∂x = E (2.32) where B models first order viscosity effects. This is basically equivalent to assuming a Voigt model for the stress-strain relation of the material. The ratio B/E is called the relaxation time. Note that both of the equations (2.31) and (2.32) have the form (jω)2 U = c2fin (ω) ∂2 U, ∂x2 (2.33) where cfin (ω) is a complex frequency-dependent wave speed given by, E for (2.31) ρ, 2 cfin (ω) = E (2.34) Bω +j , for (2.32) . ρ ρ The steady state solution to (2.33) can be found easily. Consider the characteristic equation of this differential equation, s2 + ω2 = 0, c2fin (ω) (2.35) 27 Chapter 2. Elastography Using Longitudinal Waves and denote its roots by s1 (ω) = −s2 (ω) = α(ω) + jβ(ω). Then, the steady state solution to (2.33) is of the form: U (x, jω; E, B) = A1 exp ((α(E, B) + jβ(E, B))x) +A2 exp ((−α(E, B) − jβ(E, B))x) , (2.36) where A1 and A2 are determined by boundary conditions. In our experimental setup (Fig. 2.1), the sample is confined at x = 0, i.e. u(0, t) ≡ 0, as a result U (0, jω) ≡ 0 = A1 + A2 , (2.37) and therefore A2 = −A1 . The phasor solution (2.36) in this case could be simplified to obtain, U (x, jω) = 2A1 (cos(βx) sinh(αx) +j sin(βx) cosh(αx)) , (2.38) where A1 is determined by the amplitude of the excitation at x =40mm. This solution is similar to a standing wave solution. However in a standing wave solution all the points have the same phase, whereas in equation (2.38) different spatial points have different phases. If the viscosity is negligible (α = 0) standing waves are produced. In the following subsections, the derived models will be used to propose two methods for the estimation of the viscoelastic parameters of the phantom blocks. 2.4.1 Peak to Node Method This method can be used to estimate Young’s modulus for a block of material, based on the longitudinal wave pattern observed inside the block. From (2.38) the amplitude of the solution, after some manipulation, can be found as |U (x, jω)|2 = 4A21 sinh2 (αx) + sin2 (βx) . (2.39) Now if the viscosity is small, the model (2.31) can be used and the solution to the characteristic equation (2.35) would be α + jβ = j ω E ρ , (2.40) 28 Chapter 2. Elastography Using Longitudinal Waves or, α = 0, ω β = . (2.41) (2.42) E ρ In this case, (2.39) becomes ω x . |U (x, jω)| = 2A1 sin E ρ (2.43) The distance of a node to a peak, dnp , of the wave pattern is equal to a quarter of the wavelength, therefore cfin = (4 dnp ) E = ρc2fin . ω , 2π (2.44) (2.45) In our case the node is at x =0mm and at each frequency, the peak is found by searching for a local maximum in a filtered version of the amplitude of Um (x, jω). In summary, the peak to node method estimates Young’s modulus based on the peak to node distance of a measured phasor profile Um (x, jω) using equation (2.45). But this method does not provide estimates of viscosity or relaxation time. In the next subsection another method is presented which also estimates the viscosity. 2.4.2 Model Fitting Method A more general method to characterize the viscoelastic properties of the sample block is to use the characteristic equation (2.35) to determine the complex wave speed, cfin (ω), for different excitation frequencies, ω, experimentally, and then use a parameterization of the viscoelastic properties, such as that presented in (2.34) to find the particular parameters of interest. Alternatively, one of the parameterized models (2.34) can be directly fitted to the data. We present this approach here. The viscoelastic model was chosen and the general solution (2.36) was fitted to the data by least squares optimization over the parameters A1 , A2 , E and B. Thus for each 29 Chapter 2. Elastography Using Longitudinal Waves frequency the following cost function was minimized: |Um (x, jω) − Um (x, jω; E, B)|2 . J(ω) = (2.46) x Figure 2.3 shows a typical fitted model. As can be seen the model fits the data very well. In the next section the normalized mean squared error (NMSE) is used to quantitatively ensure that the model fits the data well. 2.5 Results The proposed rheometry methods of Section 2.4 were applied to the data collected using the apparatus and methods of Section 2.3. The results for the five phantoms follow in this section. 2.5.1 Peak to Node Results The peak to node method was used to analyze the measured phasor profiles Um (x, jω) for the five phantom blocks. Note that (2.44) and (2.45) could be calculated at any frequency for which there is a node and a peak present in the phasor profile (see Fig. 2.2 (a)). As a result, this method cannot be used at very low frequencies, since no peaks occur within the 40mm height of the phantom. Figure 2.4 shows the estimated cfin (ω) for the test phantom No. 3 over the range of frequencies where this method could be used. The average speed over this frequency range, cfin , is equal to 7.5m/s and the standard deviation is equal to 0.29m/s. The Young’s modulus calculated based on this average speed is 55kPa. The same method was used to estimate longitudinal wave speed and Young’s modulus for the other four test phantoms. The results are summarized in Table 2.4. Also presented are the values obtained from conventional rheometry on the same blocks. As can be seen the results are generally in good agreement and this supports the validity of the model and thus the assumptions. 2.5.2 Model Fitting Results The model fitting approach was used to find the values of Young’s Modulus E and relaxation time for the five phantom blocks. Figure 2.5 compares the results to the conventional rheometry results. In Fig. 2.5(a) the model-fitting and rheometry values of the Young’s modulus, E, for different frequencies 30 Chapter 2. Elastography Using Longitudinal Waves (a) (b) (c) Figure 2.3: Fitted model to the measured displacement data. The data is the same as that of Fig. 2.2, which is from phantom No. 3 at 90Hz. (a) absolute value and (b) phase of the phasor profile (c) the same normalized phasor profile in the complex plane. 31 Chapter 2. Elastography Using Longitudinal Waves Figure 2.4: Estimated cfin (ω) for phantom No.3. Plotted in dark solid line is the selected region for estimation of average speed. At lower frequencies no peak occurs inside the block length. Table 2.4: Estimated values of longitudinal wave speed cfin and Young’s modulus for the test phantoms using the peak to node method peak to average node conventional wave freq. Young’s rheometry speed range modulus Young’s phantom ±std used ±std modulus±10% No (m/s) (Hz) (kPa) (kPa) error 1 4.4±0.15 30–90 19±1.5 21 11% 2 5.2±0.25 40–110 28±2.1 32 11% 3 7.5±0.29 60–150 55±4.6 56 2% 4 8.7±0.27 70–150 76±6.0 76 0% 5 12 ±0.13 90–150 140±13 127 10% 32 Chapter 2. Elastography Using Longitudinal Waves Table 2.5: Estimated values of longitudinal wave speed cfin and Young’s modulus for the test phantoms using the model fitting method model fitting conventional wave freq. Young’s rheometry speed range modulus Young’s phantom ±std used ±std modulus±10% No (m/s) (Hz) (kPa) (kPa) error 1 4.2±0.15 50–60 18±1.3 21 16% 2 5.0±0.18 40–80 26±1.9 32 18% 3 7.4±0.29 55–130 55±4.4 56 2% 4 8.7±0.34 60–135 77±6.1 76 1% 5 11.3±0.46 75–150 132±11 127 4% are plotted. The model fitting plots are limited to frequencies for which the normalized mean squared error defined by NMSE = J(ω) 2, x |Um (x, jω) − mean (Um (x, jω))| (2.47) is better than 10%. This threshold was used to ensure that the model fitted the data well enough for further conclusions to be valid. The rheometry plots are limited to 40Hz. As expected, the values of E from model fitting matches with rheometry values. The average values and standard deviations are summarized in Table 2.5. In Fig. 2.5(b) the optimized values of the relaxation times together with the values obtained from the conventional rheometry are presented. Note that both axes are in logarithmic scale in this figure. Again the data of the conventional rheometry are limited to 40Hz. As can be seen from this figure, the general trends and values of relaxation times from model fitting match that of conventional rheometry. Also a power-law decrease in the relaxation time with the increase of the frequency is observed which is consistent with the previous reports [40]. 2.6 Discussion It has been shown experimentally by Catheline et al. [36] that when a finite block of viscoelastic material is vibrated by a compressional exciter with a 33 Chapter 2. Elastography Using Longitudinal Waves (a) (b) Figure 2.5: (a) Estimated value of the Young’s modulus and (b) relaxation times for the five phantoms at different frequencies of excitation. The lower frequency graphs (4–40Hz) are from conventional rheometry and the higher frequency graphs (above 40Hz) are from model fitting. 34 Chapter 2. Elastography Using Longitudinal Waves small cross sectional area, two waves propagate in the block: a compressional wave which parts almost instantaneously at the speed of ultrasound, i.e. 1540m/s, and a longitudinally polarized slow shear wave, which follows with the shear wave speed, i.e. µ/ρ which is around 1-10m/s. The shear wave can be tracked by ultrasound speckle tracking and the displacements can be used to measure the viscoelastic properties of the block. We have shown here that when a finite block of viscoelastic material is vibrated by a compressional exciter with a large cross sectional area which covers the entire top surface of the block, instead of the longitudinally polarized slow shear wave, a slow longitudinal wave propagates with the speed cfin = E/ρ ≈ 3µ/ρ which is around 5-15m/s. This longitudinal wave can also be tracked by ultrasound speckle tracking and the displacements can be used to measure the viscoelastic properties of the block. To show how the size of the exciter affects the wave speed we conducted two different experiments with the same phantom block, namely block No. 1, similar to the experiments of [36]. In the first experiment, a plate exciter was used which covered the entire top surface of the phantom. The surface was lubricated with ultrasound gel to minimize the friction between the plate and the phantom. A single period of 50Hz sinusoid was used as a transient excitation. The exciter was an open loop voice coil. Therefore the actual mechanical vibration transferred to the phantom was not exactly a single period of 50Hz, but had some tail. One scan line (A-line) was sampled at 5kHz and the displacements tracked. The displacements at different depths are plotted as a function of time in Fig. 2.6(a). By looking at the wave front the speed is approximately 3.7m/s. This speed is 3µ/ρ which was also reported in Tables 2.4 and 2.5. The slight difference is due to a few months’ time interval between the conduction of the experiments. The properties of the phantom have shifted slightly over time. The same phantom was used for a second experiment. This time, instead of the plate, a small screw was connected to the voice coil, which acted almost like a point source. The result for this case is presented in Fig. 2.6(b). The speed from this figure is approximately 2.3m/s. This is the wave speed Catheline et al. reported in [36], i.e. the shear wave speed µ/ρ. Note that in both cases Fig. 2.6(a) and (b), just after the compression is applied, a compressional wave front is generated which propagates through the phantom at the speed of ultrasound, i.e. 1420m/s. The initial propagation of this wave front cannot be tracked by the ultrasound, since the speed of ultrasound is the same as the wave speed. This wave front probes the boundaries of the medium and gets reflected back. However the displacements caused by multiple reflections of this wave from the boundaries are 35 Chapter 2. Elastography Using Longitudinal Waves Figure 2.6: (a) The result of the excitation with a plate which covers the top surface of the phantom (b) The result of the excitation with a point source. present in the phantom. The displacements, as shown by Sandrin et al. [13], are very small and can barely be detected by our speckle tracking algorithm. Super Soft PlasticTM was shown to be a promising phantom material for elastography studies, as its Young’s (and hence shear) modulus can be adjusted over a wide range (15–150kPa), without changing its density and the speed of ultrasound cinf (Table 2.1). Figure 2.7 summarizes the Young’s modulus obtained vs. the volume portion of hardener or softener used in the fabrication process for this material. However, the relaxation time and its frequency dependent behavior is not dependent on the volume portion of Softener or Hardener, as the graphs in Fig. 2.5(b) almost overlap for the five test blocks. Therefore without any additives, this material by itself is unsuitable for creating phantoms with significant viscosity contrast. In the previous sections, measurement results for cinf (ultrasound speed) and E were presented. If a purely elastic medium is assumed, these measurements can be used to obtain all the elastic parameters. Other elastic parameters, namely Lam´e constants λ and µ, and Poisson’s ratio ν could be easily found using the relations (2.7) to (2.9). The parameters are summarized in Table 2.6 for the five phantoms. The values of the Young’s modulus from the model fitting method were used in the calculations. The speed of ultrasound cinf has little variability among different soft tissue types which are usually scanned using ultrasound imaging systems. 36 Chapter 2. Elastography Using Longitudinal Waves Figure 2.7: Obtained Young’s modulus as a function of the volume portion of Plastic Softener or Hardener added during the fabrication of the phantoms Table 2.6: Summary of the elastic properties of the five test phantoms phantom λ µ E 0.5 - ν K No (GPa) (kPa) (kPa) 1E-6× (GPa) 1 2.1 6.0 18 1.5 2.1 2 2.1 8.7 26 2.1 2.1 3 2.0 18 55 4.7 2.0 4 2.1 26 77 6.1 2.1 5 2.0 44 132 11 2.0 37 Chapter 2. Elastography Using Longitudinal Waves This fact allows the conversion of time scale of ultrasound RF-pulses to distance for ultrasound image formation. This was also true for the Super Soft PlasticTM used in our experiments. This fact has another important consequence which has been noted by other scholars as well [3]. Since soft tissue is nearly incompressible, λ µ; hence cinf is determined mainly by λ. In other words µ has no contribution. Since the variability of density is also very low (in living tissue and phantoms), λ should be almost constant. As a result, the first Lam´e constant, λ, is not an elastic parameter of interest for creating contrast in elastograms. In contrast, tissue stiffness is determined by Young’s modulus, E, which for incompressible materials is equal to 3µ, and a great variability exists here. As a result, the second Lam´e constant, µ, (or Young’s modulus, E) is the parameter of interest in elastography. From this discussion, a connection can be made between elastography methods based on shear wave excitation and methods based on longitudinal wave excitation. In shear wave excitation, shear wave speed cS = µ/ρ and hence shear modulus, µ is sought. In longitudinal excitation, Young’s modulus E = 3µ is sought. Hence these methods, although different in many aspects, essentially use the same underlying physical quantity to create elastogram contrast. Another important point is that at low frequencies of excitation the phase difference in oscillations of the different points is low because of the large wavelength. Hence it is more difficult if not impossible to estimate the viscoelastic properties by using a low frequency. We conclude that a lower limit for the frequency of excitation exists. Note that this statement only applies to the estimation of the absolute values of the viscoelastic properties using the wave equation; relative estimations are still possible though. As a matter of fact, for a static compression ∂u/∂t ≡ ∂ 2 u/∂t2 ≡ 0. In this case the wave equation (2.29) implies ∂ 2 u/∂z 2 ≡ 0. The wave equation itself will become, ρ × 0 = E × 0, (2.48) which cannot be inverted. However the strain, which is given by = ∂u/∂z is not necessarily equal to zero. Assuming a uniform stress, the strain is inversely proportional to the Young’s modulus E, and hence relative values of E can be obtained from the relative values of . If in addition, force measurement is carried out, the relative values can be turned into absolute values. It was observed that one-dimensional models would not fit the data very well at higher frequencies of excitation. It is known that a block of material acts as a wave-guide for elastic wave propagation [32], and the wave speed 38 Chapter 2. Elastography Using Longitudinal Waves depends on the dimensions of the source (the excited modes of the wave guide) as well as the frequency of excitation. Further study of these effects is needed. Finally, the assumption of a 1D model which was made for the displacements is a limiting assumption. However, the methods presented in this article are intended for rheometry purposes, i.e. characterization of the viscoelastic behavior of a test block assuming homogeneity. The traditional force-displacement rheometry methods treat the entire test block as one unit, while in the presented methods an insight into internal vibrations inside the block is obtained as well. Moreover no force measurements are carried out. The presented methods do not yield local viscoelastic properties and hence are not presented as an alternative to any of the techniques intended for that purpose. Many research groups are looking at this problem, i.e. solving the inverse problem, from different perspectives. In particular they use a more complete model, such as a 2D or 3D model, for the local inversion. This article, however, shows that longitudinal waves can also be a valid candidate for forming elastograms and extending the inverse problem methods to this type of excitation. 2.7 Summary and Conclusion In this article an experimental framework was presented for viscoelastic characterization of tissue and tissue mimicking material. Longitudinal wave excitation was used on five blocks of test phantoms and displacements were imaged by an ultrasound system. Two methods were presented for analyzing the displacements. The first method estimated the Young’s modulus based on the distance from the peak to the node of the standing wave patterns inside the block. The second method used a model fitting scheme to find the values of the Young’s modulus and the relaxation time. The results were compared to conventional rheometry results obtained using forcedisplacement measurements. For the elasticities good agreement was observed which validated the theory. For the viscosities, only the trends could be compared as the measurements were over different frequency ranges. We conclude that the speed of the longitudinal waves in finite media can be orders of magnitude lower than their speed in infinite media and in fact comparable to the speed of shear waves. It is possible to track these waves by using the ultrasound-based motion tracking techniques. The longitudinal waves cannot simply be dismissed because of their presumed high speed of propagation or their insignificant amplitude. The viscoelastic properties 39 Chapter 2. Elastography Using Longitudinal Waves of homogeneous phantom blocks can be estimated by applying the methods presented in this article on the longitudinal displacements. The methods are based on fitting a longitudinal-wave-equation-based model to the measured displacements. The estimated properties would be in agreement with the gold standard values from the conventional rheometry tests. Whether the human body or a particular tissue sample in vivo or ex vivo can be regarded as a finite medium, however, depends on the geometry of the setup and the excitation used and needs further analysis. 40 References [1] A. Sarvazyan, “Handbook of elastic properties of solids, liquids and gases,” vol. III, 2001. [2] K. Parker, L. Gao, R. Lerner, and S. Levinson, “Techniques for elastic imaging: A review,” IEEE Eng. Med. Biol. Mag., vol. 15, no. 6, pp. 52–59, 1996. [3] J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imaging elastic properties of biological tissues,” Annu. Rev. Biomed. Eng., vol. 5, pp. 57–78, April 2003. [4] K. J. Parker, L. S. Taylor, S. Gracewski, and D. J. Rubens, “A unified view of imaging the elastic properties of tissue,” J. Acoust. Soc. Am., vol. 117, no. 5, pp. 2705–2712, May 2005. [5] J. Ophir, S. K. Alam, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, C. Merritt, R. Righetti, R. Souchon, S. Srinivasan, and T. Varghese, “Elastography: Imaging the elastic properties of soft tissues with ultrasound,” Journal of Medical Ultrasonics, vol. 29, pp. 155–171, 2002. [6] R. Lerner and K. Parker, “Sonoelasticity images, ultrasonic tissue characterization and echographic imaging,” in Proceedings of the 7th European Communities Workshop, 1987, nijmegen, Netherlands. [7] R. Lerner, K. Parker, J. Holen, R. Gramiak, and R. Waag, “Sonoelasticity: Medical elasticity images derived from ultrasound signals in mechanically vibrated targets,” Acoustical Imaging, vol. 16, pp. 317– 327, 1988. [8] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: A quantitative method for imaging of elasticity of biological tissues,” Ultrason. Imaging, vol. 13, pp. 11–134, 1991. [9] O. Rudenko, A. Sarvazyan, and S. Emelianov, “Acoustic radiation force and streaming induced by focused nonlinear ultrasound in a dissipative medium,” J. Acoust. Soc. Am., vol. 99, no. 5, pp. 2791–2798, May 1996. 41 Chapter 2. Elastography Using Longitudinal Waves [10] A. Sarvazyan, O. Rudenko, S. Swanson, J. Fowlkes, and S. Emelianov, “Shear wave elasticity imaging: A new ultrasonic technology of medical diagnostics,” Ultrasound Med. Biol., vol. 24, no. 9, pp. 1419–1435, 1998. [11] L. Sandrin, M. Tanter, J.-L. Gennisson, S. Catheline, and M. Fink, “Shear elasticity probe for soft tissue with 1-d transient elastography,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 49, no. 4, pp. 436–446, April 2002. [12] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulus imaging with 2-d transient elastography,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 49, no. 2, pp. 426–435, April 2002. [13] L. Sandrin, D. Cassereau, and M. Fink, “The role of the coupling term in transient elastography,” J. Acoust. Soc. Am., vol. 115, no. 1, pp. 73–83, January 2004. [14] M. Fatemi and J. Greenleaf, “Ultrasound stimulated vibro-acoustic spectrography,” Science, vol. 280, pp. 82–85, 1998. [15] K. Hoyt, K. J. Parker, and D. J. Rubens, “Real-time shear velocity imaging using sonoelastographic techniques,” Ultrasound Med. Biol., vol. 33, no. 7, pp. 1086–1097, 2007. [16] B. Chan, J. V. Merton-Gaythrope, M. P. Kadaba, S. Zafaranloo, and D. Bryk, “Acoustic properties of polyvinyl chloride gelatin for use in ultrasonography,” Radiology, pp. 215–216, 1984. [17] T. J. Hall, M. Bilgen, M. F. Insana, and T. A. Krouskop, “Phantom materials for elastography,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 44, no. 6, pp. 1355–1365, 1997. [18] C. M. Hassan and N. A. Peppas, “Structure and applications of poly(vinyl alcohol) hydrogels produced by conventional crosslinking or by freezing/thawing methods,” Advances in Polymer Science, vol. 153, pp. 37–65, 2000. [19] E. Madsen, G. Frank, T. 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. 42 Chapter 2. Elastography Using Longitudinal Waves [20] E. L. Madsen, M. A. Hobson, G. R. Frank, H. Shi, J. Jiang, T. J. Hall, T. Varghese, M. M. Doyley, and J. B. Weaver, “Anthropomorphic breast phantoms for testing elastography systems,” Ultrasound Med. Biol., vol. 32, no. 6, pp. 857–874, 2006. [21] E. L. Madsen, M. A. Hobson, H. Shi, T. Varghese, and G. R. Frank, “Stability of heterogeneous elastography phantoms made from oil dispersions in aqueous gels,” Ultrasound Med. Biol., vol. 32, no. 2, pp. 261–270, 2006. [22] K. Zell, J. Sperl, M. Vogel, R. Niessner, and C. Haisch, “Acoustical properties of selected tissue phantom materials for ultrasound imaging,” Phys. Med. Biol., vol. 52, pp. N475–N484, 2007. [23] J. Fromageau, J.-L. Gennisson, C. Schmitt, R. L. Maurice, R. Mongrain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel mechanical properties with four ultrasound elastography methods and comparison with gold standard testings,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 54, no. 9, pp. 498–509, March 2007. [24] U. Hamhaber, F. Grieshaber, J. Nagel, and U. Klose, “Comparison of quantitative shear wave mrelastography with mechanical compression tests,” Magn. Reson. Med., vol. 49, pp. 71–77, 2003. [25] H. Mansy, J. Grahe, and R. Sandler, “Elastic properties of synthetic materials for soft tissue modeling,” Phys. Med. Biol., vol. 53, pp. 2115– 2130, 2008. [26] J. Fromageau, J.-L. Gennisson, C. Schmitt, R. L. Maurice, R. Mongrain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel mechanical properties with four ultrasound elastography methods and comparison with gold standard testing,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007. [27] M. Doyley, J. Weaver, and K. Paulsen, “Measuring the mechanical properties of soft tissue over a wide frequency range by employing the principle of time temperature super-position,” in Proceedings of the Seventh International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, Oct 2008, p. 75. [28] S. Srinivasan, T. Krouskop, and J. Ophir, “A quantitative comparison of modulus images obtained using nanoindentation with strain elastograms,” Ultrasound Med. Biol., vol. 30, no. 7, pp. 899–918, 2004. 43 Chapter 2. Elastography Using Longitudinal Waves [29] D. Kalanovic, M. P. Ottensmeyer, J. Gross, G. Buess, and S. L. Dawson, “Independent testing of soft tissue viscoelasticity using indentation and rotary shear deformations,” Stud. Health Technol. Inform., vol. 94, pp. 37–43, 2003. [30] V. Egorov, S. Tsyuryupa, S. Kanilo, M. Kogit, and A. Sarvazyan, “Soft tissue elastometer,” Med. Eng. Phys., vol. 30, pp. 206–212, 2008. [31] A. Samani, J. Bishop, C. Luginbuhl, and D. B. Plewes, “Measuring the elastic modulus of ex vivo small tissue samples,” Phys. Med. Biol., vol. 48, pp. 2183–2198, 2003. [32] H. Kolsky, Stress Waves in Solids. Inc., 1963. New York: Dover Publications, [33] K. F. Graff, Wave Motion in Elastic Solids. Oxford University Press, 1975. [34] F. Lockett, Non-linear Viscoelastic Solids. 1972. Ely House, London: London: Academic Press, [35] W. Flugge, Viscoelasticity, 2nd ed. Berlin: Springer, 1975. [36] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelkaram, and J. Culiolic, “Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: An inverse problem approach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec 2004. [37] F. A. Duck, Physical Properties of Tissue: A Comprehensive Reference Book. Academic Press, 1990. [38] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complexvalued stiffness reconstruction by for magnetic resonance elastography by algebraic inversion of the differential equation,” Magn. Reson. Med., vol. 45, pp. 299–310, 2001. [39] R. S. C. Cobbold, Foundations of Biomedical Ultrasound. University Press, 2006. Oxford [40] H. Eskandari, S. E. Salcudean, and R. Rohling, “Viscoelastic parameter estimation based on spectral analysis,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 55, no. 7, pp. 1611–1625, 2008. 44 Chapter 2. Elastography Using Longitudinal Waves [41] R. Zahiri and S. Salcudean, “Motion estimation in ultrasound images using time-domain cross-correlation with prior estimates,” IEEE Trans. Biomed. Eng., vol. 53, no. 10, pp. 1990–2000, 2006. 45 Chapter 3 Theoretical Limitations of the Elastic Wave Equation Inversion for Tissue Elastography∗ 3.1 Introduction Elastography is a fast developing field of medical imaging which strives to provide images of the mechanical properties of tissue [1–8]. The main mechanical property of interest is the elasticity or Young’s modulus E. Other mechanical properties of interest include viscosities, relaxation times, nonlinearity parameters, harmonics, poro-elasticity and Poisson’s ratio [8–22]. The majority of the elastography techniques have three components in common: • an excitation mechanism which creates quasi-static or dynamic deformations in the tissue [23–31], • a medical imaging system augmented by a motion estimation technique which is capable of providing displacement (or velocity) images of the tissue while it is being deformed [25, 32–41], and • an inversion technique which transforms the displacement images into elasticity images, the so called elastograms [25, 42–54]. For the inversion, different approaches have been devised. Most of them, however, fall under one of the two following categories: • Global inversion techniques [44–48,55–57]: These techniques typically use finite element modeling (FEM) of the medium. The goal is to ∗ A version of this chapter has been published: A. Baghani, S.E. Salcudean, R. Rohling, “Theoretical Limitations of the Elastic Wave Equation Inversion for Tissue Elastography,” J. Acoust. Soc. Am. , vol. 126, no. 3, pp. 1541-1551, Sep. 2009 46 Chapter 3. Theoretical Limitations find the best local elasticity values which would create the same displacement pattern as those observed inside the tissue, under similar boundary and excitation conditions. To find the optimum elasticity values, the forward problem is solved iteratively. In each iteration the elasticity values are adjusted until the computed displacements from the model match the measured displacement from the tissue. In addition to the measured displacements, these techniques typically require further information about the shape of the tissue, its boundary conditions, and the excitation. • Local inversion techniques [11, 25, 27–30, 37, 50–52, 58–61]: Typically in these techniques, a particular wave equation is assumed to govern the displacements. This partial differential equation (PDE) relates the spatial derivatives to the temporal derivatives of the displacement with the local elasticity as the coefficient. If the displacement measurements are available over a range of space and time, the elasticities can be found either by an algebraic inversion of the wave equation or by estimating the phase speed from the gradient of the phase. The next section is devoted to the presentation of these techniques. It is the objective of this research to: • Study the theoretical limitations of the local inversion techniques for estimating the tissue elasticity. • Study the effects of the formulation of the wave equation on the inversion results. • Study the effects of the choice of the exciter and imaging technique. To this end, the general form of the elastic wave equation in a linear continuum is studied. The hypothesis to be proven is that the phase speed, which is the final estimated quantity in many of the local inversion techniques, does not exactly represent the local mechanical properties of the tissue. The hypothesis is proven by showing that in addition to the local mechanical properties of the tissue, the phase speed depends on the geometry of the wave as well as the geometry of the tissue itself. This effect is known as the Lamb wave effect specially in the context of thin plates [62]. The practical implication is that artifacts are inevitably present in the elastograms, no matter how well the algorithm is implemented or how accurate the displacement measurements are. Nevertheless, many of the elastography techniques have been proven clinically to provide valuable information about the elasticity of the tissue. The excitation in these techniques is chosen so that the actual displacements created in the tissue satisfy the assumptions made, which in turn justifies the inversion techniques used. These meth47 Chapter 3. Theoretical Limitations ods include quasi-static constant strain compression, pulsed excitation with external exciters, and acoustic-radiation-force-based techniques. A previous study of the limitations of the local inversion techniques for estimating the tissue elasticity can be found in [63]. In that work the authors experimentally measured the displacements caused by a circular piston inside rectangular blocks of tissue mimicking material and meat specimens. In their experiments the exciter would start to vibrate at a single frequency (monochromatic excitation). Between the transient and steady state regimes of vibration, a region was identified when transient shear waves are propagating in the medium and phase speed measurements can be carried out. The deleterious effects of the medium boundaries were dealt with by rejecting the steady state regime for phase speed measurements. The effects of the boundary conditions imposed by the exciter were also studied. An approximate model, the Rayleigh–Sommerfeld’s solution, and the Green’s function were used in this analysis. They concluded that the size of the exciter affects the measured phase speed. Another conclusion was that at very low frequencies, the effects of the longitudinal wave on the phase speed cannot be neglected. The final preferred method was to use a pulsed transient excitation with a small piston to get more accurate viscoelastic measurements. The analysis of the aforementioned article is limited to the particular configuration used, which is very similar to a point source on a semi-infinite medium. In many elastography applications, to get a desired level of displacements inside tissue, it is of interest to use exciters with a larger footprint. We have shown previously that if the size of the exciter is comparable to the size of the sample, the phase speed observed is the extensional wave speed [64]. In the present article we look at the problems associated with phase speed measurements from a more general point of view, and discuss some of the phenomena that affect the accuracy of the measurements. We will show that neglecting these phenomena results in un-accounted-for artifacts in the computed elastograms. We start by a brief review of the state-of-the-art local inversion techniques in Section 3.2. Section 3.3 is a critique on the use of the terminology associated with the elastic waves. The main theoretical results of the article are presented in Section 3.4. Simulation results are presented in Section 3.5. Section 3.6 is devoted to the conclusion. 48 Chapter 3. Theoretical Limitations 3.2 Local Inversion of the Wave Equation The mechanical wave equation in an isotropic elastic continuum written in the Cartesian coordinate system is [62] ρ ∂ 2 ui ∂ = (λ + µ) ∆ + µ∇2 ui ∂t2 ∂i i = x, y, z , (3.1) where x = [x, y, z]T are the Cartesian coordinates which define the spatial location of the point, t is the time variable, u(x,t) = [ux , uy , uz ]T is the displacement field, and the coefficients ρ(x), λ(x) and µ(x) are the density, the first, and the second Lam´e constants respectively. Here, we have adopted the notation of Kolsky [62]. The dilatation ∆ is defined by ∆= ∂ux ∂uy ∂uz + + , ∂x ∂y ∂z (3.2) and determines the volumetric changes as the wave propagates. ∇2 is the Laplacian operator: ∂2 ∂2 ∂2 ∇2 = + + . (3.3) ∂x2 ∂y 2 ∂z 2 We denote the Poisson’s ratio and the Young’s modulus by ν and E respectively. The material of interest in the field of elastography is tissue in vivo. This material is found to be nearly incompressible [51, 65], hence the values of ∆ are typically very small. Other implications of the incompressibility in terms of the mechanical properties are as follows: ν ≈ 0.5 , λ µ and E ≈ 3µ . (3.4) In the rest of this article it is assumed that the relations (3.4) hold. To obtain a clinically useful image of the soft tissue, the contrast should be based on a physical quantity which has a high degree of variability among different tissue types. The density of most soft tissue is close to the density of water [66], 1000kg/m3 so ρ is not a candidate. Neither is the first Lam´e constant λ which is found to be close to the λ of water [67], i.e. 2.3GPa. On the other hand, µ and E = 3µ, which determine the stiffness of the tissue, vary over a wide range from a few to a few hundred kPa [68]. Moreover, a change in the stiffness of the tissue is often associated with pathology [69–71]. Hence elastograms are preferentially based on the tissue elasticity E = 3µ. Given the displacements as a function of time and space, the 49 Chapter 3. Theoretical Limitations most that can be recovered from equation (3.1) are the ratios λ/ρ and µ/ρ. In practice, the elastograms are based on E/ρ. Since the variability of ρ is small, the obtained contrast is almost the same as the one obtained by using the absolute value of elasticity E. A number of issues arise when medical imaging devices are used to estimate the displacement fields. In many cases it is not possible to acquire all three components of the displacement field, i.e. ux , uy , and uz . In other cases the acquired displacement field may not be available over a volume. For instance, in most cases of ultrasound elastography with 2D probes, only ux (x, y) is available and only on a single plane z =constant and not over a volume. Some customized ultrasound systems can estimate the displacements over a volume and in multiple directions [37, 41, 72]. MRI based techniques on the other hand, can acquire the full displacement field over a volume, at the expense of longer acquisition times [14, 25, 50, 51, 73–75]. Different approaches to the inversion of the wave equation have been taken. As mentioned earlier the dilatation ∆(x, t) levels in soft tissue are well beyond the accuracy of the measurement systems. Therefore the wave equation (3.1) cannot directly be inverted to obtain µ/ρ. Some state-of-theart approaches to tackle this problem are as follows: • Assuming zero dilatation or zero pressure gradient [30, 36, 37, 58, 59]: In this approach it is assumed that the dilatation is identically zero, ∆ ≡ 0, or that the pressure has no spatial variation in the medium. These assumptions reduces the wave equation to ρ ∂ 2 ui = µ∇2 ui ∂t2 i = x, y, z . (3.5) Note that for an anisotropic inversion all three components are necessary. If isotropy is assumed, measuring one component of the displacement field, for instance ux , is enough. In this case, it is natural to chose the measurement axis (the x-axis) aligned with the axis on which the wave causes the greatest displacements. If the measurement is only available on a plane, for instance ux (x, y), or on a single line, 50 Chapter 3. Theoretical Limitations ux (x), then the following assumptions are usually made ∇2 u x = ≈ ≈ ∂2 + ∂x2 ∂2 + ∂x2 ∂2 ux . ∂x2 ∂2 ∂2 + ∂y 2 ∂z 2 ∂2 ux ∂y 2 ux (3.6) • Assuming a shear wave equation [52]: In this case it is assumed that the wave is a plane shear wave. If the wave causes the greatest displacements in the x-direction, then this component of the displacement field is measured and assumed to satisfy a 1D shear wave equation, ∂ 2 ux ∂ 2 ux ρ 2 =µ 2 . (3.7) ∂t ∂y • Assuming a thin rod (extensional) wave equation [64]: When the size of a compressing exciter is comparable to the dimensions of the tissue being imaged the tissue behaves like a thin rod experiencing extensional waves. Under these conditions the component of the displacement field parallel to the direction of the excitation, say the xcomponent, satisfies ∂ 2 ux ∂ 2 ux ρ 2 =E . (3.8) ∂t ∂x2 • Writing the shear wave equation for the rotations [14, 53, 72–75]: Define the rotation w(x,t) as half the curl of the displacement field: 1 w(x,t) = ∇ × u(x, t) 2 (3.9) From (3.1) each component of the rotation field satisfies a separate shear wave equation, ρ ∂ 2 wi ∂ 2 wi = µ ∂t2 ∂i2 i = x, y, z . (3.10) Note that this approach does not make any assumptions, but it requires the measurement over a 3D volume of at least two components of the displacement field for an isotropic inversion. For an anisotropic inversion, all the components are necessary. The local inversion algorithms (with the exception of the last approach) 51 Chapter 3. Theoretical Limitations reduce the general form of the wave equation (3.1) to a 1D wave equation of the form ∂ 2 ux ∂ 2 ux ρ 2 = f (E) , (3.11) ∂t ∂x2 where f (E) is a function of the local elasticity. The quantity cph = f (E)/ρ (3.12) in this 1D equation is called the phase speed. For a unidirectional harmonic solution, ux (x, t) = a exp (jω(x/cph − t)) , (3.13) this quantity gives the propagation speed of the constant phase surfaces. Given ux (x, t), it is possible to invert equation (3.11) and find f (E)/ρ. This can be done either by directly solving for f (E)/ρ: f (E) ∂ 2 ux /∂t2 = 2 , ρ ∂ ux /∂x2 (3.14) This method is sometimes called the Algebraic Inversion of the Differential Equation (AIDE) [50]. Another method is to use the concept of the phase speed, ω f (E)/ρ = cph = ∂ , (3.15) ∂x ∠ux sometimes called the Phase Gradient method [50]. If f (E) were only a function (even an unknown function) of the local mechanical properties of the tissue, the elastograms based on the phase speed, f (E)/ρ would be free from artifacts. However as we show in Section 3.4, f (E) is not only a function of the local mechanical properties of the medium but also a function of the geometry of the wave and the geometry of the medium. An even more important result is proved; the reduction of (3.1) to (3.11) may not be possible at all since multiple values of the phase speed can coexist. 3.3 A Critique on the Terminology Associated with the Wave Equation Before we move on to our results on the phase speed we present our choice of the terminology. By the fundamental theorem of vector calculus, or Helmholtz decomposition theorem, any sufficiently smooth and decaying 52 Chapter 3. Theoretical Limitations vector field can be written as the sum of a divergence-free and a curl-free vector field [76]. The solutions of the wave equation (3.1) satisfy the conditions of the theorem, and therefore it is possible to write a vibration pattern inside an elastic medium as the sum of a divergence-free and a curl-free displacement field. • Dilatational wave: or better termed irrotational component is a component of the wave which is curl-free. This component of the wave phenomena is solely the result of voluminal changes in the medium. • Distortion wave: or better termed equivoluminal component is a component of the wave which is divergence-free. This component of the wave phenomena preserves the volumes of the infinitesimal elements. There is some confusion however, over the definition and application of the terms “shear wave”, “transverse wave”, “longitudinal wave” , and “compressional wave”. These terms are sometimes defined as follows: • Longitudinal wave: or compressional wave is a wave for which the particle velocities are parallel to the phase velocity. For a harmonic wave the particle velocity would be perpendicular to the equiphase surfaces. • Transverse wave: or shear wave is a wave for which the particle velocities are perpendicular to the phase velocity. For a harmonic wave the particle velocity would be tangent to the equiphase surfaces. These definitions are very intuitive with respect to the literal meanings of the terms. However a major problem exists with these definitions; to the best knowledge of the authors, there is no theorem proving that the quality of a wave being transversal or longitudinal (taken these definitions) is invariant with respect to the wave equation. In other words, a wave propagating in an infinite elastic medium might be purely longitudinal at one instant but not so at a later instant (even without encountering a boundary and undergoing mode conversion). This severely limits the applicability of these terms. When these definitions are taken, a connection is usually assumed between longitudinal and irrotational waves on one hand and transverse and equivoluminal waves on the other hand. The origin of this connection could be the study of the simple case of a plane wave in an infinite medium. Indeed a longitudinal plane wave propagating in an infinite medium has no equivoluminal components and consists solely of a dilatational component. This wave propagates at a phase speed of (λ + 2µ)/ρ hence the name longitudinal wave speed. Similarly a shear plane wave propagating in an infinite medium has no dilatational components and consists solely of an equivoluminal component. This wave propagates at a phase speed of µ/ρ hence 53 Chapter 3. Theoretical Limitations the name shear wave speed. A better way to define these terms, however, is to take for them the same definitions as presented for the “dilatational waves” and “distortion waves”. This is the accepted nomenclature in the physics literature, specially in the field of the electromagnetic waves [77]: • Longitudinal wave: or compressional wave is a component of the wave which is curl-free. • Transverse wave: or shear wave is a component of the wave which is divergence-free. These are the definitions we use in this article and recommend. A number of comments are in order. It is clear from the definitions that a wave can contain both the longitudinal and transversal components, or may lack one. In some cases one is not interested in the decomposition of the displacement field into these components, but rather in the study of the total displacement itself. However since each of these components separately satisfies a reduced wave equation with a definite wave speed, it is at other times useful to study them separately. The speed of the longitudinal component and the transverse component is equal to (λ + 2µ)/ρ and µ/ρ respectively. Since the wave equation can be reduced to the two aforementioned wave equations, it is evident that in an infinite elastic medium, a purely longitudinal wave will remain purely longitudinal, and a purely transverse wave will remain purely transverse. This invariance property makes the terms well-defined. The drawback is that the intuitive relationship between the directions of the particle velocity and the phase velocity is no longer valid. As a matter of fact, the quality of a wave being longitudinal or transversal could be determined solely from the displacement vector field at a single instant. The wave equation preserves this quality as the time evolves. As will become clear in the next section, the relationship between the phase speed and the longitudinal and shear wave speeds, (λ + 2µ)/ρ and µ/ρ, is a complex relationship which depends on the geometry of the wave, as well as the geometry of the medium. 54 Chapter 3. Theoretical Limitations 3.4 Phase Speed of Mechanical Waves in Elastic Mediums 3.4.1 Dependence of the Phase Speed on the Geometry of the Wave: Infinite Mediums Theorem. Consider an infinite homogeneous linear elastic medium with density ρ and Lam´e constants λ and µ. Given the angular frequency ω for a harmonic wave: • for any number c such that µ ≤c< ρ λ + 2µ , ρ (3.16) there exists a shear beam form for the harmonic wave for which the phase speed is equal to c. • for any number c such that λ + 2µ ≤ c, ρ (3.17) there exists infinitely many beam forms for the harmonic wave for which the phase speed is equal to c. These waves contain both longitudinal and shear components. The theoretical derivation which follows can be found in classical textbooks on elastic waves [62,78] as part of the solution to the wave equation in cylindrical coordinate systems. However studying the implications in terms of the phase speed, as summarized in the theorem, is novel. Proof. Consider the elastic wave equation in the cylindrical coordinate system, ∂ 2 ur ∂t2 ∂ 2 uθ ρ 2 ∂t ∂ 2 uz ρ 2 ∂t ρ ∂∆ 2µ ∂ ω ¯z ∂ω ¯θ − + 2µ , ∂r r ∂θ ∂z 1 ∂∆ ∂ω ¯r ∂ω ¯z = (λ + 2µ) − 2µ + 2µ , r ∂θ ∂z ∂r ∂∆ 2µ ∂ 2µ ∂ ω ¯r = (λ + 2µ) − (rω ¯θ ) + , ∂z r ∂r r ∂θ = (λ + 2µ) (3.18) (3.19) (3.20) where (r, θ, z) are the cylindrical coordinates and u = (ur , uθ , uz ) is the 55 Chapter 3. Theoretical Limitations displacement field. The dilatation in the cylindrical coordinates is given by, ∆= 1 ∂ 1 ∂uθ ∂uz (rur ) + + , r ∂r r ∂θ ∂z (3.21) and the rotations are given by, 2¯ ωr = 2¯ ωθ = 2¯ ωz = ∂uθ 1 ∂uz − , r ∂θ ∂z ∂ur ∂uz − , ∂z ∂r ∂ur 1 ∂ (ruθ ) − r ∂r ∂θ (3.22) (3.23) . (3.24) Consider a cylindrical wave beam propagating along the z-axis with the center of the beam on the z-axis. More specifically assume that uθ is zero, i.e. the particle displacements are confined to the rz-planes. Also assume that the wave is symmetrical around the z-axis, hence ∂/∂θ annihilates the variables. From (3.22) and (3.24) ω ¯r = ω ¯ z = 0. The reduced wave equation in this case would be, ∂ 2 ur ∂t2 ∂ 2 uz ρ 2 ∂t ρ ∂∆ ∂ω ¯θ + 2µ , ∂r ∂z ∂∆ 2µ ∂ = (λ + 2µ) − (rω ¯θ ) . ∂z r ∂r = (λ + 2µ) (3.25) (3.26) We are interested in harmonic waves propagating along the z-axis, for instance in the negative z-direction. The general form of such a wave is, ur = U (r) exp (i(kz z + ωt)) , (3.27) uz = W (r) exp (i(kz z + ωt)) , (3.28) where kz is the wave number and ω is the angular frequency. U (r) and W (r) determine the shape of the beam. Note that the phase speed of this wave is equal to cph = ω/kz . Substitution in (3.21) and (3.23) yields the expressions 56 Chapter 3. Theoretical Limitations for the dilatation and rotation: ∆(r, z, t) = ∂U (r) U (r) + + ikz W (r) exp (i(kz z + ωt)) ∂r r 2¯ ωθ (r, z, t) = ∂W (r) exp (i(kz z + ωt)) ikz U (r) − ∂r , (3.29) . (3.30) The forms (3.27) and (3.28) simplify the wave equation (3.25) and (3.26), ∂∆ + 2iµkz ω ¯θ , ∂r 2µ ∂ = ikz (λ + 2µ)∆ − (rω ¯θ ) . r ∂r −ρω 2 ur = (λ + 2µ) (3.31) −ρω 2 uz (3.32) By eliminating ω ¯ θ and ∆ between these equations, the longitudinal and transversal wave equations are derived: ∂ 2 ∆ 1 ∂∆ + + k∆ 2 ∆ = 0 , ∂r2 r ∂r ∂2ω ¯θ 1 ∂ω ¯θ ω ¯θ + − 2 + kω¯ θ 2 ω ¯θ = 0 , ∂r2 r ∂r r (3.33) (3.34) where k∆ and kω¯ θ are the geometrical beam numbers defined by k∆ 2 = kω¯ θ 2 = ω2 λ+2µ ρ 2 ω µ ρ − kz 2 , − kz 2 . (3.35) (3.36) In the case that the desired phase speed, c, satisfies (3.16) the corresponding choice of the wave number, kz = ω , c (3.37) results in a real k∆ and an imaginary kω¯ θ . In the other case (3.17) both k∆ and kω¯ θ are real. The significance of this will become clear shortly. The solution to (3.33) and (3.34) can be found by change of variables from r to k∆ r and kω¯ θ r respectively. The resulting equations are Bessel equations of the zeroth and first orders respectively. The physically mean57 Chapter 3. Theoretical Limitations ingful solutions which have bounded values at r = 0, are ∆(r, z, t) = G(z, t)J0 (k∆ r) , (3.38) ω ¯ θ (r, z, t) = H(z, t)J1 (kω¯ θ r) , (3.39) where J0 (·) and J1 (·) are Bessel functions of the first kind and of zeroth and first orders respectively. Now the forms for the dilatation and rotation are given in (3.29) and (3.30). For these forms to match the solutions (3.38) and (3.39), U (r) and W (r) should have the following forms, ∂ J0 (k∆ r) + C2 kz J1 (kω¯ θ r) ∂r = −C1 k∆ J1 (k∆ r) + C2 kz J1 (kω¯ θ r) , (3.40) i ∂ [rJ1 (kω¯ θ r)] r ∂r = C1 ikz J0 (k∆ r) + C2 ikω¯ θ J0 (kω¯ θ r) , (3.41) U (r) = C1 W (r) = C1 ikz J0 (k∆ r) + C2 where C1 and C2 are two arbitrary constants. From (3.29) and (3.30) the dilatation and rotation become ∆(r, z, t) = −2C1 (kz 2 + k∆ 2 )J0 (k∆ r) exp (i(kz z + ωt)) , (3.42) 2¯ ωθ (r, z, t) = 2iC2 (kz 2 + kω¯ θ 2 )J1 (kω¯ θ r) exp (i(kz z + ωt)) . (3.43) The wave can thus be a mixture of the longitudinal and shear components with C1 and C2 defining the respective proportions. Now if k∆ is imaginary, the Bessel functions J0 (k∆ r) and J1 (k∆ r) go to infinity as r goes to infinity. This is not physically meaningful. Thus for a phase speed c which satisfies (3.16), C1 must be zero in (3.40) and (3.41). In view of (3.42) and (3.43) this is a purely shear wave beam, which travels at the phase speed c. On the other hand if the phase speed c satisfies (3.17), both C1 and C2 can have nonzero values. Therefore infinitely many beam forms exist for which the wave travels at the phase speed c. Each of these beam forms contains both the longitudinal and shear components. This completes the proof. It follows that the phase speed depends not only on the mechanical 58 Chapter 3. Theoretical Limitations properties of the medium, but also on the geometry of the wave. Even for a purely shear wave (C1 = 0), the phase speed can have any value which is greater than or equal to µ/ρ. The shear wave speed, µ/ρ, is not the phase speed of every shear wave in an infinite medium. It is however the phase speed of the uniform plane shear waves. This issue is in addition to the main drawback of the use of the phase speed [63] for tissue characterization; namely that the phase speed can not be defined when multiple waves are traveling in different directions, for instance when there are reflections of the wave from the boundaries. 3.4.2 Dependence of the Phase Speed on the Geometry of the Medium: Wave Guides In the previous section the medium was assumed to be infinite in size. This assumption is also made in many of the elastography techniques [63, 79]. However no part of the human body is infinite in size. In this section we present some classical results on the wave guides and study their implications in the field of elastography. The wave guides are infinite in at least one direction and finite in at least another. Therefore they cannot model the tissue behavior accurately either. However they can be considered as an intermediate step in moving from the analysis of an infinite medium to a bounded medium. As such the insight gained from studying them is useful in understanding and designing elastography systems. The wave guide we study is an infinitely long cylindrical rod. Choose the cylindrical coordinate system with the axis of the cylinder on the z-axis. As in the previous section, we are interested in harmonic waves propagating along the z-axis, i.e. along the axis of the cylinder, which are symmetrical around the z-axis. The same analysis applies and the wave pattern given by (3.27) and (3.28) with U (r) and W (r) given by (3.40) and (3.41) satisfies the wave equation with phase speed c, provided that kz is chosen to satisfy ω/kz = c. However in this case the boundary conditions impose restrictions on the permissible values of c. The boundary of the cylinder is free from stresses. The expressions for 59 Chapter 3. Theoretical Limitations stresses in the cylindrical coordinate system are, ∂ur , ∂r ∂ uθ 1 ∂ur +r = µ r ∂θ ∂r r ∂ur ∂uz = µ . + ∂z ∂r σrr = λ∆ + 2µ σrθ σrz (3.44) , (3.45) (3.46) By the assumptions σrθ vanishes everywhere. The boundary conditions require that σrr and σrz vanish on the surface of the cylinder. Denote the radius of the cylinder by a. Substitution of the expressions for ur and uz from (3.27), (3.28), (3.40), and (3.41) translates the boundary conditions into the following two equations, C1 2µ ∂2 ∂r2 2C1 kz ∂ ∂r λρω 2 J0 (k∆ a) λ + 2µ ∂ +2C2 µkz J1 (kω¯ θ r) = 0 , ∂r r=a J0 (k∆ r) − r=a (3.47) J0 (k∆ r) r=a +C2 2kz 2 − ρω 2 µ J1 (kω¯ θ a) = 0 . (3.48) Note that these equations depend only on the ratio of C1 /C2 . Eliminating this variable between the two, a single equation is obtained for the unknown kz . This equation is called the Pochhammer frequency equation [62]. For a given material ρ, µ, λ and a are known. So the equation basically describes the relationship between the wave number kz and angular frequency ω of the wave, or equivalently the relationship between the phase speed c and ω. As it turns out, unlike the infinite medium, not every phase speed is possible for a given ω. The permissible speeds must satisfy the Pochhammer frequency equation. Substitution into either of the equations (3.47) or (3.48) determines the ratio C1 /C2 , i.e. the proportions of the longitudinal and transversal waves that should be mixed together to obtain that particular phase speed. 60 Chapter 3. Theoretical Limitations Table 3.1: Mechanical property values used for the simulation ρ λ µ kg/m3 GPa kPa 1000 2.3 10 3.5 3.5.1 Simulation Results Shear Waves in a Tissue Mimicking Infinite Medium To gain some insight into the theoretical results which were presented in the previous section, we simulate the shear wave beam forms in an infinite medium for different angular frequencies and phase speeds. For a chosen pair of angular frequency ω and phase speed c, the wave number kz is found from (3.37). We choose the mechanical properties of the medium to match those found in the soft living tissue such as the breast. These values are listed in Table 3.1. Given the mechanical properties of the medium, k∆ and kω¯ θ are found from (3.35) and (3.36) respectively. In the next step U (r) and W (r) are found from (3.40) and (3.41), by setting C1 = 0 and choosing an arbitrary value for C2 . Note that C1 is chosen to be zero to obtain a purely shear wave and the value of C2 only determines the overall amplitude and phase of the wave, but does not change the waveform. If the time evolution of the wave is desired, U (r) and W (r) can be substituted in (3.27) and (3.28) to obtain ur (r, z, t) and uz (r, z, t). √ As shear waves can propagate at any phase speed above µ/ρ = 10m/s, the choice of the phase speed, c, becomes arbitrary. We simulated the beam forms using two chosen values of c = 5m/s and c = 12m/s. The results are shown in Fig. 3.1 and 3.2 respectively. The first and the second column in these figures show the displacement components ur ((r, θ, z), t) and uz ((r, θ, z), t) respectively at t = 0. The three rows correspond to increasing frequencies of the wave 40Hz, 65Hz, and 100Hz. Note that because of the uniqueness of the solution to the wave equation, if a hypothetical planar exciter could be built which would create the same harmonic motion as ur ((r, θ, z), t) and uz ((r, θ, z), t) over an infinite cross section of the medium, for instance at z = 0, the wave forms due to this cylindrically symmetric infinite exciter, in the steady state, would be the same as those depicted in these figures. The figures show how changing the excitation pattern of such an exciter can result in a completely different phase speed, even if the frequency of the excitation is not changed. 61 Chapter 3. Theoretical Limitations Figure 3.1: Infinite medium; shear beam patterns for three different frequencies all sharing the same phase speed of 5m/s. (a), (b), (c) radial components of the displacement (d), (e), (f) z-component of the displacement at 40Hz, 65Hz, and 100Hz respectively. Note that only a cubic portion of the medium near the axis of symmetry is shown; The medium extends in all directions to infinity. 62 Chapter 3. Theoretical Limitations Figure 3.2: Infinite medium; shear beam patterns for three different frequencies all sharing the same phase speed of 12m/s. (a), (b), (c) radial components of the displacement (d), (e), (f) z-component of the displacement at 40Hz, 65Hz, and 100Hz respectively. Note that only a cubic portion of the medium near the axis of symmetry is shown; The medium extends in all directions to infinity. 63 Chapter 3. Theoretical Limitations Table 3.2: Mechanical property values used for the simulation ρ λ µ a kg/m3 GPa kPa m 1000 2.3 10 0.05 3.5.2 Waves in a Tissue Mimicking Cylindrical Rod The Pochhammer equation has been studied for metallic rods such as steel beams. To gain some insight into this equation when dealing with the living tissue, we consider a simple example here. We choose the mechanical properties of the medium to match those found in the living tissue. These values are listed in Table 3.2. Using these values, the Pochhammer equation was solved numerically using MATLAB for different frequencies. Figure 3.3 shows the plots of the values of c obtained for each ω. At higher frequencies, the equation has multiple roots (modes), thus the multiples plots in this figure. Mode 0 has a constant phase speed c = µ/ρ, i.e. the shear wave speed for all frequencies. However substitution of the shear wave speed into (3.40) and (3.41) results in U (r) = W (r) = 0. Therefore this mode is a trivial solution. As a matter of fact assuming a phase speed equal to the shear wave speed results in vanishing displacements, independent of the material properties. The implication is that no (axis symmetric) wave can travel along the cylindrical rod with the shear wave speed. Figure 3.3: The first four modes of the cylindrical bar 64 Chapter 3. Theoretical Limitations At low frequencies (below 38Hz) only mode 1 is present (see Fig. 3.3). Therefore the low frequency waves can only travel at a single speed. This is in contrast to the case of the infinite medium studied before, in which the waves, independent of their frequency, could travel at a range of speeds. The low frequency speed is equal to E/ρ which is the familiar value obtained from the thin rod approximation theory [64]. This also shows the theoretical justification behind the assumptions made in (3.8). As the frequency goes higher, other modes start to appear. At high frequencies (above 38Hz) multiple waves of the same frequency can propagate simultaneously, each with a different phase speed. For instance at a frequency of 100Hz, three phase speeds of 3.0339m/s, 3.719m/s, and 5.630m/s are possible. The corresponding beam patterns are shown in Fig. 3.4. In this case it may not be possible to recover a phase speed from studying the displacement patterns inside the material. To see this more clearly, assume that a harmonic exciter vibrating at a frequency of 100Hz is placed at infinity on the rod and has caused the following wave pattern: ur (t) = 3Uc1 (r) exp (i(207.1z + 2π100t)) +5Uc2 (r) exp (i(168.94z + 2π100t)) +8Uc3 (r) exp (i(111.60z + 2π100t)) , (3.49) uz (t) = 3Wc1 (r) exp (i(207.1z + 2π100t)) +5Wc2 (r) exp (i(168.94z + 2π100t)) +8Wc3 (r) exp (i(111.60z + 2π100t)) , (3.50) where the Uci and Wci are the solutions presented in Fig. 3.4 for the three phase speeds. Because of the linearity of the wave equation, any linear combination of the solutions presented in Fig. 3.4 is a solution. In particular the above linear combination with the coefficients 3, 5, and 8 caused by the particular exciter used, is a solution. It is not hard to verify by substitution that neither the direct inversion method (3.14), nor the phase gradient method (3.15) results in a meaningful mechanical property for the homogeneous cylinder. The implication in the context of elastography is that at higher frequencies, multiple modes appear in the measured displacements which makes it impossible to recover a single phase speed and determine the mechanical properties based on that. 65 Chapter 3. Theoretical Limitations Figure 3.4: Permissible beam patterns (modes) in an infinite cylindrical rod for waves having a frequency of 100Hz. Three (and only three) phase speeds are possible. (a), (b), (c) radial components of the displacement (d), (e), (f) z-component of the displacement for phase speeds of 3.0339m/s, 3.719m/s, and 5.630m/s respectively. 66 Chapter 3. Theoretical Limitations 3.6 Conclusion The propagation of elastic waves trapped in finite media, such as soft tissue, is a complex phenomenon even without considering the viscous and nonlinear effects. The term wave is generally applied to any vibration pattern inside the tissue. However the phase velocity cannot be defined for all the cases of the vibration patterns. Usually multiple waves traveling in different directions superimpose to create the vibration pattern. In this case the phase velocity might be well-defined for each wave, but not for the resultant of the superposition. Even in the simple case where a single frequency wave is traveling in a single direction, the geometry of the wave and the medium create multiple permissible phase speeds, and thus make it difficult to recover a meaningful phase speed for the traveling wave. The parameters used in our simulations were chosen to be close to those of the living tissue. The diameter of the cylindrical rod was chosen to be 10cm which is in the same size range as human organs. Yet multiple phase speeds and modal behavior start to appear at frequencies as low as 38Hz. This frequency is in the range of frequencies used in many dynamic excitation elastography techniques. One remedy seems to be the use of the natural decomposition of the vibrations into the longitudinal (dilatational) and transverse (shear) components. Each of these components always satisfies its own wave equation with its own wave speed which is independent of the geometry (note that this is not the phase speed however). Taking the divergence of the measured displacement field yields the dilatation which satisfies the dilatational wave equation with wave speed (λ + 2µ)/ρ, ρ ∂2∆ = (λ + 2µ)∇2 ∆ . ∂t2 (3.51) Since living tissue is incompressible, the dilatations will be very minute and impossible to detect by taking the spatial derivatives of the measured displacements from any of the currently available imaging techniques. Even if accurate measurements were possible, inverting this wave equation would result in the local values of the longitudinal wave speed, (λ + 2µ)/ρ, being known. However the change of this parameter is quite small in the soft tissue (from 1400m/s for fat to 1540m/s for muscle). Therefore the contrast of the formed elastogram would be minimal. On the other hand, taking the curl of the displacement field results in 67 Chapter 3. Theoretical Limitations the rotation fields, ∇ × (ux , uy , uz )T = (¯ ωx , ω ¯y , ω ¯ z )T each of which satisfies a shear wave equation with wave speed ρ ∂2ω ¯i = µ∇2 ω ¯i ∂t2 i = x, y, z . (3.52) µ/ρ, (3.53) Since living tissue is incompressible, E ≈ 3µ and inverting these equations results in local information for tissue elasticity E. Moreover since these equations are naturally 1D, it is possible to use the phase speed in their context. From this discussion it is concluded that the only theoretically flawless method of imaging tissue elasticity is by taking the curl of the displacement fields. This requires the measurement of the displacements in all three directions over the volume of interest. In cases where this is not feasible, such as ultrasound elastography, artifacts will always be present. These artifacts are not just due to noisy measurements, poor algorithms, or imperfect setups but due to the theoretical limitations of the systems themselves. There are some measures which could be taken to minimize the artifacts. The main goal should be to create a single wave propagating in one direction at a single speed: • The excitation amplitude should be chosen small enough to reduce the reflected wave amplitudes to the level of other measurement noises. The damping present in the tissue helps by reducing the amplitude of the reflected waves from the boundaries of the tissue and bones. • The frequency of excitation should be kept to a minimum to prevent the appearance of higher modes. However lower frequencies result in lower damping, and therefore a trade off exists here. • The excitation pattern should be chosen so as to excite mainly the lowest mode. From an engineering point of view, however, many factors beyond the analysis presented here affect the success of an elastography system. Many of the devised elastography techniques are proven to provide clinically valuable images and some of them have even found their way into the market. Quasistatic constant strain, pulsed excitation with external exciters, acousticradiation-force-based and other techniques are each designed to make the displacements created in the clinical setting match the assumptions made about them. This results in the validity of their inversion techniques and 68 Chapter 3. Theoretical Limitations their valuable elastograms. The presence of the artifacts, by itself, does not mean that an elastography system should be dismissed. After all, each of the available medical imaging modalities has its own artifacts, and yet they are very well proven to be useful in diagnosis and treatment. 69 References [1] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: A quantitative method for imaging of elasticity of biological tissues,” Ultrasonic Imaging, vol. 13, pp. 111–134, 1991. [2] K. Parker, L. Gao, R. Lerner, and S. Levinson, “Techniques for elastic imaging: A review,” IEEE Engineering in Medicine and Biology Magazine, vol. 15, no. 6, pp. 52–59, 1996. [3] L. Gao, K. J. Parker, R. M. Lerner, and S. Levinson, “Imaging of the elastic properties of tissue-a review,” Ultrasound Med. Biol., vol. 22, no. 8, pp. 957–977, 1996. [4] J. Ophir, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, R. Righetti, and T. Varghese, “Elastographic imaging,” Ultrasound Med. Biol., vol. 26, pp. S23–S29, 2000. [5] A. Sarvazyan, Handbook of Elastic Properties of Solids, Liquids and Gases. Academic Press, 2001, vol. III, ch. Elastic Properties of Soft Tissues. [6] J. Ophir, S. K. Alam, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, C. Merritt, R. Righetti, R. Souchon, S. Srinivasan, and T. Varghese, “Elastography: Imaging the elastic properties of soft tissues with ultrasound,” Journal of Medical Ultrasonics, vol. 29, pp. 155–171, 2002. [7] J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imaging elastic properties of biological tissues,” Annual Review Biomedical Engineering, vol. 5, pp. 57–78, April 2003. [8] K. J. Parker, L. S. Taylor, S. Gracewski, and D. J. Rubens, “A unified view of imaging the elastic properties of tissue,” J. Acoust. Soc. Am., vol. 117, no. 5, pp. 2705–2712, May 2005. [9] T. Varghese, J. Ophir, and T. Krouskop, “Nonlinear stress-strain relationships in tissue and their effect on the contrast-to-noise ratio in elastograms,” Ultrasound Med. Biol., vol. 26, no. 5, pp. 839–851, 2000. 70 Chapter 3. Theoretical Limitations [10] E. E. Konofagou, T. P. Harrigan, J. Ophir, and T. A. Krouskop, “Poroelastography: Imaging the poroelastic properties of tissues,” Ultrasound Med. Biol., vol. 27, no. 10, pp. 1387–1397, 2001. [11] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity and viscosity from measurement of shear wave speed dispersion,” J. Acoust. Soc. Am., vol. 115, no. 6, pp. 2781–2785, June 2004. [12] J. Bercoff, M. Tanter, M. Muller, and M. FinkJ, “The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 51, no. 11, pp. 1523–1536, 2004. [13] R. Righetti, J. Ophir, S. Srinivasan, and T. A. Krouskop, “The feasibility of using elastography for imaging the poisson’s ratio in porous media,” Ultrasound Med. Biol., vol. 30, no. 2, pp. 215–228, 2004. [14] R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink, “Viscoelastic shear properties of in vivo breast lesions measured by mr elastography,” Magnetic Resonance Imaging, vol. 23, pp. 159–165, 2005. [15] R. Righetti, J. Ophir, and T. A. Krouskop, “A method for generating permeability elastograms and poisson’s ratio time-constant elastograms,” Ultrasound Med. Biol., vol. 31, no. 6, pp. 803–816, 2005. [16] S. Chen, R. Kinnick, J. F. Greenleaf, and M. Fatemi, “Difference frequency and its harmonic emitted by microbubbles under dual frequency excitation,” Ultrasonics, vol. 44, pp. 123–126, 2006. [17] R. Sinkus, J. Bercoff, M. Tanter, J.-L. Gennisson, C. E. Khoury, V. Servois, A. Tardivon, and M. Fink, “Nonlinear viscoelastic properties of tissue assessed by ultrasound,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 53, no. 11, pp. 2009–2018, 2006. [18] S. Chen, R. R. Kinnick, J. F. Greenleaf, and M. Fatemi, “Harmonic vibro-acoustography,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 54, no. 7, pp. 1346–1351, July 2007. [19] X. Jacob, S. Catheline, J.-L. Gennisson, C. Barri´ere, D. Royer, and M. Fink, “Nonlinear shear wave interaction in soft solids,” J. Acoust. Soc. Am., vol. 122, no. 4, pp. 1917–1926, October 2007. 71 Chapter 3. Theoretical Limitations [20] A. Thitaikumar, T. A. Krouskop, B. S. Garra, and J. Ophir, “Visualization of bonding at an inclusion boundary using axial-shear strain elastography: a feasibility study,” Phys. Med. Biol., vol. 52, pp. 2615– 2633, 2007. [21] J. Gennisson, M. R´enier, S. Catheline, C. Barri´ere, J. Bercoff, M. Tanter, and M. Fink, “Acoustoelasticity in soft solids: Assessment of the nonlinear shear modulus with the acoustic radiation force,” J. Acoust. Soc. Am., vol. 122, no. 6, pp. 3211–3219, December 2007. [22] A. Thitaikumar, L. M. Mobbs, C. M. Kraemer-Chant, B. S. Garra, and J. Ophir, “Breast tumor classification using axial shear strain elastography: a feasibility study,” Phys. Med. Biol., vol. 53, pp. 4809–4823, 2008. [23] R. Lerner and K. Parker, “Sonoelasticity images, ultrasonic tissue characterization and echographic imaging,” in Proceedings of the 7th European Communities Workshop, 1987, nijmegen, Netherlands. [24] R. Lerner, K. Parker, J. Holen, R. Gramiak, and R. Waag, “Sonoelasticity: Medical elasticity images derived from ultrasound signals in mechanically vibrated targets,” Acoustical Imaging, vol. 16, pp. 317– 327, 1988. [25] R. Sinkus, J. Lorenzen, D. S. amd M. Lorenzen, M. Dargatz, and D. Holz, “High-resolution tensor mr elastography for breast tumor detection,” Physics Med. Biol., vol. 45, pp. 1649–1664, 2000. [26] R. Souchon, L. Soualmi, M. Bertrand, J.-Y. Chapelon, F. Kallel, and J. Ophir, “Ultrasonic elastography using sector scan imaging and a radial compression,” Ultrasonics, vol. 40, pp. 867–871, 2002. [27] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulus imaging with 2-d transient elastography,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 49, no. 4, pp. 426–435, April 2002. [28] L. Sandrin, M. Tanter, J.-L. Gennisson, S. Catheline, and M. Fink, “Shear elasticity probe for soft tissues with 1-d transient elastography,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 49, no. 4, pp. 436– 446, April 2002. [29] J.-L. Gennisson, S. Catheline, S. Chaffa¨ı, and M. Fink, “Transient elastography in anisotropic medium: Application to the measurement of 72 Chapter 3. Theoretical Limitations slow and fast shear wave speeds in muscles,” J. Acoust. Soc. Am., vol. 114, no. 1, pp. 536–541, July 2003. [30] J. Bercoff, R. Sinkus, M. Tanter, and M. Fink, “Supersonic shear imaging: A new technique for soft tissue elasticity mapping,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 51, no. 4, pp. 396–409, 2004. [31] M. W. Urban and J. F. Greenleaf, “Harmonic pulsed excitation and motion detection of a vibrating reflective target,” J. Acoust. Soc. Am., vol. 123, no. 1, pp. 519–533, Jan 2008. [32] F. Yeung, S. F. Levinson, D. Fu, and K. J. Parker, “Feature-adaptive motion tracking of ultrasound image sequences using a deformable mesh,” IEEE Trans. Med. Imaging, vol. 17, no. 6, pp. 945–956, 1998. [33] E. Konofagou, T. Varghese, J. Ophir, and S. Alam, “Power spectral strain estimators in elastography,” Ultrasound Med. Biol., vol. 25, no. 7, pp. 1115–1129, 1999. [34] E. E. Konofagou and J. Ophir, “Precision estimation and imaging of normal and shear components of the 3d strain tensor in elastography,” Phys. Med. Biol., vol. 45, pp. 1553–1563, 2000. [35] T. Varghese, E. Konofagou, J. Ophir, S. Alam, and M. Bilgen, “Direct strain estimation in elastography using spectral cross-correlation,” Ultrasound Med. Biol., vol. 26, no. 9, pp. 1525–1537, 2000. [36] M. Fink, L. Sandrin, M. Tanter, S. Catheline, S. Chaffai, J. Bercoff, and J. Gennisson, “Ultra high speed imaging of elasticity,” in IEEE Ultrasonics Symposium, 2002, pp. 1811–1820. [37] J. Bercoff, R. Sinkus, M. Tanter, and M. FinkJ, “3d ultrasound-based dynamic and transient elastography: First in vitro results,” in IEEE Ultrasonics Symposium, 2004, pp. 28–31. [38] K. Hoyt, F. Forsberg, and J. Ophir, “Investigation of parametric spectral estimation techniques for elasticity imaging,” Ultrasound Med. Biol., vol. 31, no. 8, pp. 1109–1121, 2005. [39] ——, “Comparison of shift estimation strategies in spectral elastography,” Ultrasonics, vol. 44, pp. 99–108, 2006. [40] ——, “Analysis of a hybrid spectral strain estimation technique in elastography,” Phys. Med. Biol. 51, vol. 51, pp. 197–209, 2006. 73 Chapter 3. Theoretical Limitations [41] J.-L. Gennisson, T. Deffieux, R. Sinkus, P. Annic, M. Pernot, F. Cudeiro, G. Montaldo, M. Tanter, M. Fink, and J. Bercoff, “A 3d elastography system based on the concept of ultrasound-computed tomography for in vivo breast examination,” in IEEE Ultrasonics Symposium, 2006, pp. 1037–1040. [42] A. J. Romano, J. J. Shirron, and J. A. Bucaro, “On the noninvasive determination of material parameters from a knowledge of elastic displacements: Theory and numerical simulation,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 45, no. 3, pp. 751–759, May 1998. [43] C. Sumi, A. Suzuki, and K. Nakayama, “Estimation of shear modulus distribution in soft tissue from strain distribution,” IEEE Trans. Biomed. Eng., vol. 42, no. 2, pp. 193–202, February 1995. [44] K. Raghavan and A. E. Yagle, “Forward and inverse problems in elasticity imaging of soft tissues,” IEEE Transactions on Nuclear Science, vol. 41, no. 4, pp. 1639–1648, August 1994. [45] F. Kallel and M. Bertrand, “Tissue elasticity reconstruction using linear perturbation method,” IEEE Trans. Med. Imaging, vol. 15, no. 3, pp. 299–313, June 1996. [46] M. Doyley, P. Meaney, and J. Bamber, “Evaluation of an iterative reconstruction method for quantitative elastography,” Phys. Med. Biol., vol. 45, pp. 1521–1540, 2000. [47] P. E. Barbone and J. C. Bamber, “Quantitative elasticity imaging: What can and cannot be inferred from strain images,” Phys. Med. Biol., vol. 47, pp. 2147–2164, 2002. [48] A. A. Oberai, N. H. Gokhale, and G. R. Feij´oo, “Solution of inverse problems in elasticity imaging using the adjoint method,” Inverse Problems, vol. 19, pp. 297–313, 2003. [49] D. Fu, S. Levinson, S. Gracewski, and K. Parker, “Non-invasive quantitative reconstruction of tissue elasticity using an iterative forward approach,” Phys. Med. Biol. 45 (2000) 14951509, vol. 45, pp. 1495–1509, 2000. [50] A. Manduca, T. Oliphant, M. Dresner, J. Mahowald, S. Kruse, E. Amromin, J. Felmlee, J. Greenleaf, and R. Ehman, “Magnetic resonance elastography: Non-invasive mapping of tissue elasticity,” Medical Image Analysis, vol. 5, pp. 237–2540, 2001. 74 Chapter 3. Theoretical Limitations [51] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complexvalued stiffness reconstruction by for magnetic resonance elastography by algebraic inversion of the differential equation,” Magnetic Resonance in Medicine, vol. 45, pp. 299–310, 2001. [52] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelkaram, and J. Culiolic, “Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: An inverse problem approach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec 2004. [53] B. Robert, R. Sinkus, B. Larrat, M. Tanter, and M. Fink, “A new rheological model based on fractional derivatives for biological tissues,” in IEEE Ultrasonics Symposium, 2006, pp. 1033–1036. [54] X. Zhang and J. F. Greenleaf, “Estimation of tissues elasticity with surface wave speed (l),” J. Acoust. Soc. Am., vol. 122, no. 5, pp. 2522– 2525, November 2007. [55] A. A. Oberai, N. H. Gokhale, M. M. Doyley, and J. C. Bamber, “Evaluation of the adjoint equation based algorithm for elasticity imaging,” Phys. Med. Biol., vol. 49, pp. 2955–2974, 2004. [56] M. M. Doyley, S. Srinivasan, S. A. Pendergrass, Z. Wu, and J. Ophir, “Comparative evaluation of strain-based and model-based modulus elastography,” Ultrasound Med. Biol., vol. 31, no. 6, pp. 787–802, 2005. [57] M. M. Doyley, S. Srinivasan, E. Dimidenko, N. Soni, and J. Ophir, “Enhancing the performance of model-based elastography by incorporating additional a priori information in the modulus image reconstruction process,” Phys. Med. Biol., vol. 51, pp. 95–112, 2006. [58] J. Bercoff, S. Chaffai, M. Tanter, L. Sandrin, S. Catheline, M. Fink, J. L. Gennisson, and M. Meunier, “In vivo breast tumor detection using transient elastography,” Magnetic Resonance in Medicine, vol. 29, no. 10, pp. 1387–1396, 2003. [59] J. Bercoff, M. Muller, M. Tanter, and M. Fink, “Study of viscous and elastic properties of soft tissue using supersonic shear imaging,” in IEEE Ultrasonics Symposium, 2003, pp. 925–928. [60] J. Fromageau, J.-L. Gennisson, C. Schmitt, R. L. Maurice, R. Mongrain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel mechanical properties with four ultrasound elastography methods and 75 Chapter 3. Theoretical Limitations comparison with gold standard testing,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007. [61] A. Romano, J. Bucaro, P. Abraham, and S. Dey, “Inversion methods for the detection and localization of inclusions in structures utilizing dynamic surface displacements,” in Proceedings of SPIE Vol. 5503, vol. 5503, 2004, pp. 367–374. [62] H. Kolsky, Stress Waves in Solids. Inc., 1963. New York: Dover Publications, [63] S. Catheline, F. Wu, and M. Fink, “A solution to diffraction biases in sonoelasticity: The acoustic impulse technique,” J. Acoust. Soc. Am., vol. 105, no. 5, May 1999. [64] A. Baghani, H. Eskandari, T. Salcudean, and R. Rohling, “Measurement of tissue elasticity using longitudinal waves,” in Proceedings of the 7th International conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, October 2008, p. 103. [65] F. A. Duck, Physical Properties of Tissue: A Comprehensive Reference Book. Academic Press, 1990. [66] M. M. Burlew, E. L. Madsen, J. A. Zagzebski, R. A. Bahjavic, and S. W. Sum, “A new ultrasound tissue-equivalent material,” Radiology, vol. 134, pp. 517–520, 1980. [67] S. A. Goss, R. L. Johnston, and F. Dunn, “Comprehensive compilation of empirical ultrasonic properties of mammalian tissues,” J. Acoust. Soc. Am., vol. 64, no. 2, pp. 423–457, August 1978. [68] A. Sarvazyan, O. Rudenko, S. Swanson, J. Fowlkes, and S. Emelianov, “Shear wave elasticity imaging: A new ultrasonic technology of medical diagnostics,” Ultrasound Med. Biol., vol. 24, no. 9, pp. 1419–1435, 1998. [69] M. Walz, J. Teubner, and M. Georgi, “Elasticity of benign and malignant breast lesions,” in Imaging, Application and Results in Clinical and General Practice, 1993, p. 56, eight International Congress on the Ultrasonic Examination of the Breast. [70] W.-C. Yeh, P.-C. Li, Y.-M. J. amd Hey-Chi Hsu, P.-L. Kuo, M.-L. Li, P.-M. Yang, and P. H. Lee, “Elastic modulus measurements of human liver and correlation with pathology,” Ultrasound Med. Biol., vol. 28, no. 4, pp. 467–474, 2002. 76 Chapter 3. Theoretical Limitations [71] T. Krouskop, T. Wheeler, F. Kallel, B. Garra, and T. Hall, “Elastic moduli of breast and prostate tissues under compression,” Ultrasonic Imaging, vol. 20, no. 4, pp. 260–274, October 1998. [72] M. Muller, J.-L. Gennisson, T. Deffieux, R. Sinkus, P. Annic, G. Montaldo, M. Tanter, and M. Fink, “Full 3d inversion of the viscoelasticity wave propagation problem for 3d ultrasound elastography in breast cancer diagnosis,” in IEEE Ultrasonics Symposium, 2007, pp. 672–675. [73] L. Huwart, F. Peeters, R. Sinkus, L. Annet, N. Salameh, L. C. ter Beek, Y. Horsmans, and B. E. V. Beers, “Liver fibrosis: Non-invasive assessment with mr elastography,” NMR in Biomedicine, vol. 19, pp. 173–179, March 2006. [74] R. Sinkus, M. Tanter, S. Catheline, J. Lorenzen, C. Kuhl, E. Sondermann, and M. Fink, “Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography,” Magnetic Resonance in Medicine, vol. 53, no. 2, pp. 372–387, 2005. [75] 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, pp. 1135– 1144, 2007. [76] R. Baierlein, “Representing a vector field: Helmholtz’s theorem derived from a fourier identity,” Am. J. Phys., vol. 63, pp. 180–182, 1995. [77] F. Rohrlich, “Causality, the Coulomb field, and Newton’s law of gravitation,” Am. J. Phys., vol. 70, no. 4, pp. 411–414, April 2002. [78] K. F. Graff, Wave Motion in Elastic Solids. Oxford University Press, 1975. Ely House, London: [79] L. Sandrin, D. Cassereau, and M. Fink, “The role of the coupling term in transient elastography,” J. Acoust. Soc. Am., vol. 115, no. 1, pp. 73–83, January 2004. 77 Chapter 4 The Influence of the Boundary Conditions on Longitudinal Wave Propagation in a Viscoelastic Medium∗ 4.1 Introduction The propagation of an acoustic wave in soft tissue is governed by local mechanical properties. For tissue characterization, such as in the field of elastography, the goal is to estimate the mechanical properties from measurements of the wave propagation. In such an inverse problem, the unknown mechanical parameters can be identified by analyzing the dynamic motion of the medium. Of the two types of waves that can propagate in a viscoelastic environment, shear waves have received significantly more attention than longitudinal waves in the field of tissue characterization. For a shear wave, particle displacements are in a direction perpendicular to the direction of the wave propagation, while for a longitudinal wave, the particle motions are in the same direction as the wave propagation. Shear waves are usually produced inside tissue either by external transverse motion of the boundary [4, 29], or through remotely induced acoustic radiation force impulse (ARFI) [2, 10, 20, 24]. The resulting shear wave can be tracked by ∗ A version of this chapter has been published: H. Eskandari, A. Baghani, S. E. Salcudean and R. Rohling, “The Influence of the Boundary Conditions on Longitudinal Wave Propagation in a Viscoelastic Medium,” Phys. Med. Bio., vol. 54, no. 13, pp. 3997-4017, Jun. 2009 78 Chapter 4. Influence of Boundary Conditions ultrasound or by magnetic resonance elastography (MRE) and the shear modulus can be reconstructed from the velocity of the wave front [19]. Longitudinal waves are assumed to have a very high propagation speed which makes them untraceable by ultrasound or other relatively low framerate data acquisition methods. Therefore, the effect of longitudinal wave is usually eliminated by applying the curl operator to the displacement field [25]. We have recently shown that it is possible to track the motion of the medium induced by longitudinal waves [1, 8]. Observations of wave propagation in tissue-mimicking phantoms indicate that a longitudinal wave can travel in finite media at a much lower speed than the speed of ultrasound, which makes it possible to track the wave front by means of ultrasound. Using the ultrasound imaging modality and a longitudinal excitation scheme, the relaxation behavior of phantoms and breast tissue were observed and quantified in [12] and [13]. In [22], a method for generating low-speed longitudinal waves was proposed. A low-frequency compressional excitation was applied by two rods that were placed on both sides of an ultrasound probe. As a result, a wave could be seen in front of the probe, having the same speed as the shear wave. This phenomenon has been explained in [22] and [23] as the effect of superposition of the shear wave fronts of the two exciters, superimposing in the middle line and producing a longitudinally polarized shear wave through mode conversion. Based on a Green’s function analysis in [21], this effect was associated with the presence of a coupling term which gives rise to a longitudinal wave propagating at the shear wave speed. In this paper we analyze the effect of the boundary conditions on the propagation speed of longitudinal waves in a finite medium. It will be shown here by exploiting the three-dimensional (3D) wave equation that a longitudinal wave may travel at a very high speed such as ultrasound or at speeds as low as that of the shear wave. In fact, for nearly incompressible soft materials with appropriate boundary conditions, a thin rod situation can be approximated where a longitudinal wave travels at a speed of nearly 1.7 times that of shear wave. The theoretical results will be verified by harmonic and transient 3D finite element simulations and the effect of viscosity will be explored. A 3D linear dynamic FEM model is adopted as a reference framework for the analysis and to validate the analytical results. Transient and harmonic simulations are performed. Experiments on a homogeneous incompressible soft material with different boundary conditions validate the longitudinal wave speeds predicted by the analytical derivations and simulations. As a result of this study, besides the shear wave, tracking the longitudinal wave in ultrasound viscoelasticity imaging can be used as a 79 Chapter 4. Influence of Boundary Conditions new method to investigate the viscoelastic behavior of soft tissue. 4.2 Governing Equations in a Linear Viscoelastic Medium In the Cartesian coordinate system, the 3D displacement vector u = [ux , uy , uz ]T and the strain tensor = [ xx , yy , zz , γxy , γxz , γyz ]T have the following relationship: ∂x 0 0 0 ∂y 0 0 0 ∂z u . = (4.1) ∂y ∂x 0 ∂z 0 ∂x 0 ∂z ∂y In a linear isotropic viscoelastic medium, the stress tensor σ is related to the time-dependent strain tensor through a combination of the Lam´e constants λ and µ, and viscosity constants λ and µ : λ + 2µ λ λ 0 0 0 λ λ + 2µ λ 0 0 0 λ λ λ + 2µ 0 0 0 σ = 0 0 0 µ 0 0 0 0 0 0 µ 0 + 0 0 0 0 0 µ λ + 2µ λ λ 0 0 0 λ λ + 2µ λ 0 0 0 λ λ λ + 2µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ ∂t , (4.2) where ∂t is the time-derivative. Here, µ and µ are often referred to as shear elasticity and shear viscosity, and λ is known as the elongational or compressional viscosity. The equation of motion for a body that is subjected to an external force 80 Chapter 4. Influence of Boundary Conditions field f , is: ∂x σxx + ∂y σxy + ∂z σxz ρ ∂t2 u = ∂x σxy + ∂y σyy + ∂z σyz + f , ∂x σxz + ∂y σyz + ∂z σzz (4.3) where u is the displacement vector as a function of time and position, and ρ is the density. By inserting from (4.1) and (4.2) into (4.3), a hyperbolic partial differential equation (PDE) results which describes the wave propagation in a linear viscoelastic medium [16]: ρ ∂t2 u = (λ + µ)∇∇ · u + µ∇2 u + (λ + µ )∇∇ · (∂t u) + µ ∇2 (∂t u) + f . (4.4) The gradient and divergence operators are denoted by (∇) and (∇·), respectively. ∇2 is the vector Laplacian operator. Note that in the Cartesian coordinate system ∇·u = xx + yy + zz , which is a measure of the volumetric change in the medium and is called the dilatation. The Young’s modulus E is defined as: E = 2(1 + ν)µ , (4.5) λ . 2(λ + µ) (4.6) where, ν= ν is the Poisson’s ratio, which takes values between 0 and 0.5. A value close to 0.5 denotes incompressibility of the material, as is the case for most soft tissues. Different assumptions have been made in the literature to reduce the number of viscosity parameters. For example, if the viscous effects are assumed to follow the elastic effects in the material, then λ /λ = µ /µ. However, since there has been little justification for such similar behavior of viscosity and elasticity, this assumption seems implausible [26]. Another assumption as made by [2] and [25] is that the compressional viscosity λ is negligible compared to the shear viscosity. To our knowledge, there is no experimental validation for this assumption for incompressible soft tissue. Another approach is to translate the idea of elastic incompressibility into the viscosity domain by assuming that the bulk viscosity is negligible, so that the fluid incompressibility condition holds. When the divergence of flow is zero, the fluid is considered to be incompressible and the term associated with the bulk viscosity in the Navier-Stokes equations will vanish [18]. In this case, the power dissipation in the medium is only a diviatoric dissipation 81 Chapter 4. Influence of Boundary Conditions associated with the distortion or shape-change, rather than dilatational dissipation that involves a change of volume [18]. The bulk viscosity is defined as κ = λ + 2µ /3 which can be set to zero to obtain: 2 λ =− µ . 3 (4.7) Bulk viscosity is known to be negligible compared to shear viscosity for rubber-like materials as well as some liquids and gases [16]. We proceed with the assumption that both conditions of elastic incompressibility (ν ≈ 0.5 in (4.6)) and fluid incompressibility (4.7) are satisfied in soft tissue. In the next section, the 3D viscoelastic wave equation is simplified for different sets of boundary conditions. The effective speed of the longitudinal wave propagation will then be extracted from the resulting representation of the wave equation. 4.3 Wave Equation and the Propagation Speed In a purely elastic infinite medium, the wave equation (4.4) becomes: ρ ∂t2 u = (λ + µ)∇∇ · u + µ∇2 u + f . (4.8) The general solution of equation (4.8) can be decomposed into a divergencefree and a curl-free component [16]. The divergence-free part tends to maintain the volume and is referred to as the shear or s-wave which travels at a speed of: µ cs = . (4.9) ρ The curl-free component has no rotation within infinitesimal blocks of material and is called the dilatational or p-wave which travels at a higher speed: cp = λ + 2µ . ρ (4.10) While plane shear or dilatational waves in an infinite medium travel at the designated speeds mentioned above, the influence from the boundary conditions may give rise to waves that travel at other velocities. A few examples are presented in the following sections. 82 Chapter 4. Influence of Boundary Conditions 4.3.1 Longitudinal Wave in a Finite Medium Approximating a Thin Rod Let us first consider the case where the excitation is mainly compressional in the x direction. If the exciter is large relative to the cross-section such that a plane-wave is produced in the medium, and the sides of the region have free boundaries, all the stress components, except for σxx , are negligible. Note that for the viscoelastic case, the strain distribution will be affected by the viscous stresses as well as the elastic stresses, however due to the compressional excitation and the plane-wave assumption, σxx is dominant and the approximation should be valid. Using this condition in (4.2), (λ + µ) + (λ + µ )∂t σxx = µ + µ ∂t (3λ + 2µ) + (3λ + 2µ)∂t xx . (4.11) For the interior region that is not experiencing any internal force, (4.3) is reduced to ∂x σxx = ρ ∂t2 ux . This can be substituted into (4.11), together with (4.1), to obtain the viscoelastic wave equation in the x -direction: ρ (λ + µ) + (λ + µ )∂t ∂t2 ux = µ + µ ∂t (3λ + 2µ) + (3λ + 2µ)∂t ∂x2 ux . (4.12) This is a third order PDE with respect to the variable t. The case of longitudinal wave propagation in a viscoelastic thin rod has been addressed more than seven decades ago by Thompson [26]. The conditions described in this section in which non-axial stress components become zero approximate a thin rod situation. Therefore, following [26], by taking the Laplace transform in the x domain (L{f (x)} = F (s)) and the Fourier transform in the t domain (F{f (t)} = fˆ(jω)), the characteristic equation of (4.12) becomes: −ω 2 ρ + s2 (µ + jωµ ) [(3λ + 2µ) + jω(3λ + 2µ )] =0. [(λ + µ) + jω(λ + µ )] (4.13) There are three roots associated with this characteristic equation with respect to ω. However, considering that at relatively low frequencies (λ+µ) (λ + µ )ω for incompressible materials, the root that corresponds to the denominator is very large and may be ignored [26]. Also, with the assumption of zero bulk viscosity, 3λ + 2µ = 0. Therefore, equation (4.13) reduces to: −ω 2 ρ + s2 (µ + jωµ )(3λ + 2µ) =0. (λ + µ) (4.14) 83 Chapter 4. Influence of Boundary Conditions Given the incompressibility of the material and the definition of Young’s modulus, this can be written as −ω 2 ρ + s2 (E + jωη) = 0, which results in the following PDE: (4.15) ρ∂t2 ux = E∂x2 ux + η∂x2 ∂t ux . η is called the coefficient of normal viscosity, where, η ≈ 3µ . (4.16) The PDE in equation (4.15) can be discretized using the well-known distributed system of several mass-spring-damper (MSD) elements which are connected to each other in series [27]. Based on (4.16), in order for the 1D model to follow the behavior of the 3D model in its best sense, the equivalent viscosity of the 1D model should be three times the shear viscosity. This will be taken into account in the following simulations where behaviors of different models are compared. Also, given that the Young’s modulus in nearly incompressible materials is approximately three times the shear modulus, the ratio between the observed viscous and elastic effects, entitled relaxation-time, will be independent of the modality of the excitation, whether it is shear or longitudinal. Equation (4.15) indicates that a longitudinal wave is present in the medium, propagating at a much lower speed than cp . In the purely elastic case, without the influence of viscosity, the speed is: c0 = E . ρ (4.17) This value, being approximately 1.7 times the shear wave speed for incompressible materials, denotes that a low-speed longitudinal wave can be induced in the medium by applying certain boundary conditions. Note that if the viscosity is also taken into account, the effective modulus from (4.15) will be E 2 + ω 2 η 2 instead of E, which results in a speed equal to c ≈ c0 (1 + ω 2 τ 2 /4), where τ = η/E is the relaxation-time. However, with the small viscosity and relaxation-time in soft tissue, this correction term is negligible at low frequencies and equation (4.17) can be used to calculate the wave speed. 4.3.2 Longitudinal Wave Induced by a Point Source Exciter Suppose a point source excitation is applied axially to a finite sample of homogeneous material in which radial symmetry can be assumed. In this 84 Chapter 4. Influence of Boundary Conditions fx r ur è uè ux x Figure 4.1: Displacement vector in the cylindrical coordinate system. case, it is more convenient to analyze the problem in cylindrical coordinates as shown in Fig. 4.1. The stress-strain relationship in the cylindrical coordinate system is the same as (4.2). Similar to the analysis in Section 4.2, in cylindrical coordinates where u = [ux , ur , uθ ]T and = [ xx , rr , θθ , γxr , γxθ , γrθ ]T , the strain tensor is defined as follows: = ∂x 0 0 ∂r 1 r ∂θ 0 0 ∂r 1 r ∂x 1 r 0 ∂θ 0 0 1 r ∂θ 0 ∂x (∂r − 1r ) u . (4.18) The equation of motion in the x -direction in this case becomes: 1 1 ρ ∂t2 ux = ∂x σxx + ∂r σxr + σxr + ∂θ σxθ + fx . r r (4.19) Note that if all of the shear stress components are zero, the resulting wave equation would be that of the p-wave in a viscoelastic medium. For the interior, where fx = 0, (4.19) can be expanded as: ρ ∂t2 ux = [(λ + 2µ) + (λ + 2µ )∂t ] ∂x ( −2(µ + µ ∂t )∂x ( +(µ + µ ∂t )∂r xr + 1 r (µ + µ ∂t ) xx rr + xr + + rr θθ ) 1 r (µ + θθ ) + µ ∂t )∂r xθ . (4.20) 85 Chapter 4. Influence of Boundary Conditions It can be shown that equation (4.20) can be simplified to: ρ ∂t2 ux = [(λ + 2µ) + (λ + 2µ )∂t ] ∂x ( −2(µ + µ ∂t ) 1r xx + rr + θθ ) [∂r (rω ¯ θ ) − ∂θ (¯ ωr )] , (4.21) where ω ¯ r and ω ¯ θ are rotations about the r and θ axes: 2¯ ωθ = ∂x ur − ∂r ux , and, 2¯ ωr = 1 ∂θ ux − ∂x uθ . r (4.22) To obtain a relationship between the radial and angular strains on the axis of excitation, l’Hˆopital’s rule can be utilized at r → 0. According to this theorem, if two functions are differentiable in a neighborhood of a given point while both of them converge to zero, the indeterminate limit of their quotients would be equal to the quotient of their derivatives at that point. By applying this rule to the angular strain, one gets: θθ |r=0 1 = lim ur = lim ∂r ur = r→0 r r→0 rr |r=0 . (4.23) Also, due to radial symmetry in this case, σθθ = 0, which according to (4.2), leads to: (λ + 2µ) + (λ + 2µ )∂t θθ + (λ + λ ∂t )( xx + rr ) =0. (4.24) Inserting from (4.23) into (4.24), at the axis of excitation: rr = θθ = −(λ + λ ∂t ) 2 [(λ + µ) + (λ + µ )∂t ] xx . (4.25) Inserting (4.25) into (4.21), noting that ω ¯ r is zero on the excitation axis and assuming that ω ¯ θ and ∂r ω ˆ θ are also zero at r → 0, the following results: ρ (λ + µ) + (λ + µ )∂t ∂t2 ux = µ + µ ∂t (λ + 2µ) + (λ + 2µ)∂t ∂x2 ux . (4.26) With the same approach as in section 4.3.1 and given the large value of λ compared to µ, µ and λ in incompressible materials, the wave equation simplifies to: ρ∂t2 ux = µ∂x2 ux + µ ∂x2 ∂t ux . (4.27) Therefore, given radial symmetry in the medium, a point-source exciter 86 Chapter 4. Influence of Boundary Conditions will generate a longitudinal wave that travels at the shear speed: µ , ρ c0 = (4.28) The assumption that ω ¯ θ and ∂r ω ˆ θ vanish at r → 0 is verified in Section 4.5 by observing the same wave as predicted above in simulations and experiments. Also, note that using (4.25), the axial stress will be: σxx = [(λ + 2µ) + (λ + 2µ )∂t ] = xx + (λ + λ ∂t )( rr (µ+µ ∂t )[(3λ+2µ)+(3λ +2µ )∂t ] xx [(λ+µ)+(λ +µ )∂t ] ≈E xx + η ∂t xx + θθ ) , (4.29) where η is as in equation (4.16). In other words, the apparent elastic constant in the axial direction is the Young’s modulus rather than the shear modulus, however the velocity of the wave is determined by the shear modulus according to (4.28). In the following section, it will be shown that for the case of a homogeneous block of material, the displacement spectra and particularly the resonance frequency of the system can be used to estimate the wave speed in the medium. The effect of viscosity on the resonance frequency and the speed will be formulated. The analysis in the following section will justify a method that will be used in Section 4.5 in order to estimate the wave speed from the displacement spectrum. 4.4 Determining the Wave Speed from the Spectra of the Displacements In this section, the effect of the viscosity on the wave equation, and in particular on the displacement spectrum, will be analyzed. It will be shown that the wave speed can be estimated from the spectrum of the displacements in a viscoelastic medium. The result from this section will be used to estimate the wave speed from the simulation results in Section 4.5, where the effect of the boundary conditions on the wave speed is explored. As shown in [1], if the medium is bounded from one end, the solution to (4.15) in the frequency domain is in the following form: u ˆ(x, ω) = A0 [cos(βx)sinh(αx) + jsin(βx)cosh(αx)] , (4.30) 87 Chapter 4. Influence of Boundary Conditions where u ˆ is the Fourier transform of ux , A0 is the amplitude of the vibration at frequency ω and α and β can be calculated by taking the Laplace transform of (4.15) in the x domain and the Fourier transform in the t domain [1]: (E + jωη)s2 + ω 2 ρ = 0 . Defining relaxation-time by τ , in soft tissue, ωτ = ωη/E the solution to (4.31) can be written as: (4.31) 1 for ω < 100Hz, jω jω jωτ s = ±(α + jβ) = ± √ ≈ ± (1 − ), c0 2 c0 1 + jωτ (4.32) where c0 is defined in (4.17). Hence, ω2τ , 2c0 ω β≈ . c0 α≈ (4.33) If the excitation is applied at x = L with a magnitude of one at all frequencies, it can be shown through simple trigonometric and hyperbolic transformations, that the magnitude of the displacement in (4.30) is: |ˆ u(x, ω)| = sinh2 (αx) + sin2 (βx) . sinh2 (αL) + sin2 (βL) (4.34) Equation (4.34) is a 2D function of ω and x, from which the wave speed c0 can be obtained. Two points on this functions are of particular interest to us: (i) The frequency f1 at which the displacement at x = L/2 is larger compared to other locations. (ii) The frequency f2 at which the displacement at x = L/2 has its maximum value compared to other frequencies. Note that for a purely elastic medium, these two points will be identical. In this way, the wave speed can be estimated by finding either the frequency that maximizes the displacement at the middle point or the frequency for which the middle point has the highest displacement magnitude compared to other locations. Figure 4.2 shows the displacement magnitude in equation (4.34) as a function of frequency and depth for a low viscosity and a high viscosity medium. A Young’s modulus of 10kP a and density of 1000kg/m3 were assumed, thus c0 = 3.16m/s and f0 = c0 /(2L) = 26.4Hz. First, with 88 Chapter 4. Influence of Boundary Conditions a viscosity of 5P as, α and β were calculated from (4.32). The displacement magnitude and the equal contour plots are shown in Figs. 4.2(a) and 4.2(b) where the frequencies as noted by case (i) and case (ii) are shown as f1 and f2 respectively. It can be seen that if the viscosity is very low, the two cases are almost the same and both f1 and f2 occur very close to f0 . In the second case, a viscosity of 30P as was assumed and the displacement magnitudes were calculated. The magitude and contour plots are depicted in Figs. 4.2(c) and 4.2(d). In this case, f1 = 33.2Hz and f2 = 24.4Hz, so using f2 to reconstruct the wave speed would be more accurate. In the following, we will explore how the wave speed can be reconstructed from case (i) using f1 and from case (ii) using f2 . The errors that each approximation produces will be quantified. The more accurate method will be adopted for wave speed analysis in the remaining sections of the paper. 4.4.1 Case (i) According to (i), ω1 should be found such that ∂x |ˆ u(x, ω1 )| = 0 at x = L/2. By differentiating (4.34) with respect to x, ∂x |ˆ u(x, ω)| = αsinh(2αx) + βsin(2βx) . sinh2 (αL) + sin2 (βL) (4.35) This equation should be zero at [x, ω] = [L/2, ω1 ]. Therefore, 1 ω1 τ sinh(αL) + sin(βL) = 0 . 2 (4.36) Based on the assumption that ωτ 1 and that ωL c0 , we have sinh(αL) ≈ αL. As a consequence, it can be shown from (4.36) that sin(βL) ≈ 0 while βL is not close to zero, so βL ≈ π. Thus, sin(βL) ≈ π − βL. From (4.36) and (4.33): ω13 τ 2 L ω1 L +π− =0. (4.37) 4c0 c0 If ω1 = 2πf1 , the wave speed can be estimated as follows: c0 = 2f1 L 1 − π 2 f12 τ 2 . (4.38) Apparently in this case, if one uses 2f1 L instead of (4.38), the wave speed will be over-estimated by a small fraction. 89 Chapter 4. Influence of Boundary Conditions 0 1 50 1 2 40 2 3 30 4 20 5 10 6 10 20 30 40 Frequency (Hz) 50 Depth (cm) Depth (cm) 0 3 4 5 6 60 20 (a) f1, f2 40 60 Frequency (Hz) (b) 0 0 1 1 2 1 3 4 Depth (cm) Depth (cm) 1.5 2 3 4 0.5 5 6 5 10 20 30 40 Frequency (Hz) (c) 50 60 6 20 f2 f1 40 60 Frequency (Hz) (d) Figure 4.2: The magnitude of the displacement in a 1D model of the wave equation with the purely elastic wave speed of 3.16m/s, length of 6cm, and Young’s modulus of 10kP a. (a) and (b) are the magnitude and contour plots for the case in which the viscosity is 5P as. The frequencies denoted by case (i) and case (ii) are almost equal. (c) and (d) show the the plots for the case that the viscosity of the medium is 30P as. In this case, f1 over-estimates the actual value, while f2 results in a relatively small under-estimation. 4.4.2 Case (ii) In this case, ω2 should be found such that ∂ω |ˆ u(x, ω)| = 0 at [x, ω] = [L/2, ω2 ]. By taking the derivative of (4.34) with respect to ω and setting it to zero for x = L/2, the following will be obtained: (ω2 τ sinh(αL) + sin(βL)) sinh2 (αL) + sin2 (βL) = αL βL 2 (ω2 τ sinh(2αL) + sin(2βL)) sinh2 ( ) + sin2 ( ) . 2 2 (4.39) 90 Chapter 4. Influence of Boundary Conditions Using the characteristics of the trigonometric and hyperbolic functions, it can be shown that: sinh2 (αL) + sin2 (βL) = (cosh(αL) + cos(βL)) (cosh(αL) − cos(βL)) , (4.40) and, αL βL cosh(αL) − cos(βL) sinh2 ( ) + sin2 ( )= . (4.41) 2 2 2 Inserting (4.40) and (4.41) into (4.39) gives: (ω2 τ sinh(αL) + sin(βL)) (cosh(αL) + cos(βL)) = ω2 τ sinh(2αL) + sin(2βL) , (4.42) which can be simplified to: sin(βL) = ω2 τ sinh(αL) . (4.43) Using the same assumptions as noted in Case (i), the wave speed can be calculated as: c0 = 2f2 L 1 + π 2 f22 τ 2 . (4.44) where ω2 = 2πf2 . If a small amount of under-estimation is tolerable, the wave speed can be approximated as 2f2 L. By comparing (4.38) and (4.44), it can be seen that f2 < f1 , therefore the second-order term in (4.44) is smaller than the second-order term in (4.38). As a result, it would be a more accurate estimation to approximate the wave speed as 2f2 L rather than 2f1 L, if the relaxation-time of the medium is unknown. Also, note that although equation (4.15) has been used here to analyze the accuracy of the two cases above, the same argument is valid for other forms of the wave equation, for example when shear modulus replaces Young’s modulus as in (4.27). 4.5 4.5.1 Results Dynamic Finite Element Analysis Numerical analysis of the three-dimensional (3D) wave equation is possible through discretization of equation (4.4) using the Ritz-Galerkin method [6]. The resulting finite element model can be used to explain the static and dynamic deformations in the medium in a discretized platform. An equation can be obtained to describe the relationship between the nodal displace91 Chapter 4. Influence of Boundary Conditions Excitation X Z 6 cm Y T 4 cm [0, 0, 0] 4 cm Figure 4.3: The geometry of the 3D region used for FEM simulations. ments, boundary conditions and local material properties: Ku(t) + B u(t) ˙ + Mu ¨(t) = f (t) , (4.45) where K, B and M are the stiffness, viscosity and mass matrices of the elements, u(t) is the nodal displacement vector and f (t) is a vector containing the boundary conditions. This equation can be written in the frequency domain as follows: (K + jωB − ω 2 M )ˆ u = fˆ . (4.46) A method to generate an accurate viscosity matrix based on a linear viscoelastic model of the continuum has been proposed in [9]. An implicit solution for the transient case of (4.45) with step-size ∆t can be obtained as [3]: ¯ Ku(t + ∆t) = f¯(t) , ¯ =K+ B + M , K ∆t ∆t2 2M B M f¯(t) = f (t) + u(t) + ( − )u(t − ∆t) . ∆t2 2∆t ∆t2 (4.47) A finite element simulation can be used to verify the theory presented in Section 4.3. In particular, numerical tests can facilitate analysis of the propagation of the fast and slow longitudinal waves in soft materials. For this purpose, a region of size 6 × 4 × 4(cm3 ) was meshed by eight-node brick elements, having a total of 25 × 7 × 7 nodes. The bottom surface at x = 0 was confined from axial movement in the x direction and the excitation 92 Chapter 4. Influence of Boundary Conditions was applied axially at the top surface, as shown in Fig. 4.3. Although the derivations in section 4.3 have been based on either radial symmetry or thin rod assumptions, the Cartesian-based rectangular FEM model utilized here does not satisfy those conditions. However, the geometry used here can provide additional insight into the longitudinal wave propagation in a non-ideal situation. In this manuscript, the middle line refers to the line in the phantom with {y, z} = [2, 2]cm, the middle point or middle node refers to the point or node on top of the phantom with {x, y, z} = [6, 2, 2]cm, and the center point is the point at {x, y, z} = [3, 2, 2]cm. 4.5.2 Simulation of the Wave Propagation The 3D region in Fig. 4.3 was assumed to consist of a homogeneous material with a Young’s modulus of 10kP a, a viscosity of 10P as and a density of 1000kg/m3 . A Poisson’s ratio of ν = 0.495 is often considered for finite element analysis of soft tissues [11]. With the given Young’s modulus and Poisson’s ratio, it can be shown that λ = 334.1kP a and µ = 3.3kP a. Three settings for the boundary conditions were simulated. The longitudinal wave patterns were compared to the theoretical results in Section 4.3 for the following cases, when: (a) The excitation was applied to all of the nodes at the top surface of the phantom. The four vertical sides of the phantom were free, thus bulging was allowed. (b) The excitation was applied to all of the nodes at the top surface of the phantom. The four vertical sides of the phantom were confined with slip boundaries such that bulging was not allowed. (c) The excitation was applied to a single node at the center of the top surface, i.e. the middle node. The four vertical sides of the phantom were free, thus bulging was allowed. In all these cases, a slip boundary condition was imposed on the bottom of the phantom, with the phantom being axially confined. In case (a), an elastic longitudinal wave is expected to be generated in the medium, traveling at a speed designated by (4.17). In case (b), bulging is not possible due to the boundary conditions. Therefore, a dilatational wave will be generated at a speed determined by (4.10). Finally, in case (c), a point source compressional exciter is modeled within the resolution of the finite element mesh. As shown in section 4.3.2, the speed of the longitudinal component of the generated wave in this case is approximated by (4.28). These analytical derivations will be tested in the following simulations. 93 Chapter 4. Influence of Boundary Conditions TF Magnitude at the Center Point TF Phase at the Center Point 1 10 0 0 ∠ H(f) [rad] |H(f)| 10 −1 10 −2 10 Free sides Bounded sides Point source −2 −4 Free sides Bounded sides Point source −6 −3 10 1 2 10 10 Frequency (Hz) (a) 1 2 10 10 Frequency (Hz) (b) Figure 4.4: The magnitude (a) and the phase (b) of the transfer function between the center point of the phantom and the middle point at the top. The solid line shows the case where the four vertical sides of the phantom are free to move and the dashed lines are for the case where the sides of the phantom are confined from bulging. The gray dash-dot line shows the transfer function for a one-node excitation at the middle point. For each case, a frequency sweep analysis was performed by solving equation (4.46) over a range of frequencies. The transfer function between the displacements at the center point and the middle point was then measured by dividing the displacements spectra at the two locations. As explained in section 4.4, by determining the frequency at which the peak of the magnitude of this transfer function occurs, the wave speed can be estimated given the small relaxation-time of the medium. The transfer function magnitudes for the three cases above are shown in Fig. 4.4(a) and the phases are plotted in Fig. 4.4(b). The peak of the three transfer functions occur at 24.4Hz, 153Hz and 14.6Hz for cases (a), (b) and (c), respectively. These frequencies can be multiplied by a wavelength of 12cm, which is twice the axial extent of the phantom, to obtain the longitudinal wave speed. Therefore, the estimated wave speed for these cases would be respectively, 2.93m/s, 18.36m/s and 1.75m/s. The actual values as computed from equations (4.10), (4.17) and (4.28) are 3.16m/s, 18.38m/s and 1.82m/s. Under different loading conditions, the effective viscosity of the medium would change. In case (a), the coefficient of normal viscosity can be calculated from (4.16) to be 30P as, which should be divided by the Young’s modulus to get a relaxation-time of τa = 1ms. In case (b), the effective viscous modulus would be λ +2µ , which is a counterpart of the effective elastic 94 Chapter 4. Influence of Boundary Conditions Effect of Viscosity on the Magnitude of TF 1 10 |H(f)| 0 10 1 Pas 5 Pas 10 Pas 20 Pas 100 Pas −1 10 1 10 2 10 Frequency (Hz) Figure 4.5: Effect of the viscosity on the magnitude of the transfer function between the center point and the middle point. The excitation was applied to the entire top surface while the sides were free to bulge. modulus λ + 2µ. Given a bulk viscosity of zero, λ + 2µ = 4µ /3 = 13.3P as. Thus, the relaxation-time in this case would be τb = 39µs. Note that in reality where λ is in the order of 109 , the viscosity for p-wave propagation would be even less significant. Lastly, for case (c), the viscosity would be equal to µ according to (4.28), which can be divided by the shear modulus to get a relaxation-time of τc = 1ms. Having the relaxation-time of the medium, the estimation of the wave speed can be enhanced by using equation (4.44). The resulting approximations of the longitudinal wave speed are, 2.95m/s, 18.37m/s and 1.76m/s for the three cases, respectively. The presence of viscosity has a crucial effect on the peak value of the transfer function magnitude. Higher viscosity results in higher dissipation and less power transfer in the medium. Figure 4.5 shows the effect of changing the viscosity on the magnitude of the transfer function for the same model as explained above. The elastic wave was generated in this case by compressing the top of the phantom and allowing the sides to move freely as in case (a). Note that the location of the peak and thus the speed of propagation do not change significantly, while the attenuation depends on the damping value. To show the transient behavior of the model, a step displacement with 1mm amplitude was applied to the top of the simulated phantom. Figure 4.6 shows how the step propagates along the middle line of the phantom. The previous three cases were explored here as well. First, an elastic wave was generated in a medium with free sides by applying the excitation to the top surface. Based on the peak of the displacement waveform in Fig. 4.6(a), it 95 Chapter 4. Influence of Boundary Conditions takes approximately 21ms for the step to travel 6cm and reach the bottom. The speed is therefore estimated as 2.9m/s. Also from Fig. 4.6(a), the resonance of the phantom in response to the step is at 24.4Hz which complies with the previous result. The magnified version of this figure is shown in Fig. 4.6(b) for a few nodes on the middle line. This figure clearly illustrates the way that information travels in the medium. It can be seen that a fast traveling wave propagates in the phantom, which reaches the bottom after 2.8ms. This corresponds to a speed of 21.4m/s which is comparable to the p-wave speed in the medium, i.e. 18.38m/s. It should be noted that determining the temporal location of the peak only gives an estimate of the wave-front rather than rendering it accurately. For the second test as in case (b), a compressional step excitation was applied to the top surface while the material was bounded from the sides, as shown in Fig. 4.6(c) and 4.6(d). Since the p-wave is the only wave that is produced in the medium, a magnified version in Fig. 4.6(d) only shows the presence of only one wavefront at a speed of 18.2m/s. Finally, the middle node on the top of the phantom was excited by a step function. The middle line displacements are depicted in Fig. 4.6(e). While the elastic wave speed is estimated at 2.1m/s, a p-wave can be identified in the magnified version in Fig. 4.6(f) with a speed of 21.4m/s. Note that in all cases, the resonances that are produced in the phantom are consistent with the peaks of the transfer functions in Fig 4.4(a). Table 4.1 gives the theoretical as well as the estimated wave speeds. The estimated values that were calculated based on the peak of the transfer function for the center node, before and after viscosity compensation are shown. Also, the wave speed that is recovered from the temporal peak of the wave-front are reported for the three cases of boundary conditions. 4.5.3 Experiments A soft homogeneous PVC (polyvinyl chloride) phantom was used to experimentally validate the speeds that were derived in Section 4.3. The shape of the phantom was as depicted in Fig. 4.3, but the dimensions were 4 × 3 × 3(cm3 ) where the axial extent was 4cm. Independent rheometry measurements resulted in a Young’s modulus of 17kP a and a density of 1000kg/m3 as reported in [1]. This will result in an elastic longitudinal wave speed of 4.1m/s and a shear wave speed of 2.4m/s. An ultrasound linear probe and a Sonix RP machine (Ultrasonix Corp., Richmond, BC, Canada) were used to capture RF data from a single A-line in the middle of the phantom at 5000 frames per second in the Doppler mode. The sides 96 Chapter 4. Influence of Boundary Conditions of the phantom were free to bulge, the transducer was placed at the bottom and a step excitation was applied from the top. The average strain in the phantom due to the step was approximately 0.1%. The axial displacements were estimated from the RF data using the peak search algorithm [7]. In the first experiment, the exciter was a plate entirely covering the top surface of the block. The displacements of different nodes along the acquisition line are shown in Fig. 4.7(a) where every displacement waveform is shifted up to the initial position of the block to which it corresponds. By following the peak of the wavefront in time, an estimate of the speed can be obtained. The peak location of the wavefront is plotted versus time in Fig. 4.7(b) for the first half of the phantom from the top. Fitting a line to the data reveals a speed of 4.0m/s for the wave travelling longitudinally in this case. One can also estimate the speed from the resonance frequency of the phantom. As seen in Fig. 4.7(a) the period of the vibration is 19.0ms which denotes a resonance frequency of 52.6Hz. Given that this frequency corresponds to a wavelength twice as large as the axial length of the phantom, the wave speed can be estimated as c = 2 × 4cm × 52.6Hz = 4.2m/s. In the second experiment, a point source exciter was approximated by using a small screw with a cross-section of 2mm2 to vibrate the phantom. The displacement waveforms are shown in Fig. 4.7(c), where the initial front Excitation (a) Flat excitation Free sides (b) Flat excitation Bounded (c) Point source Theory (m/s) Estimation (m/s) Spectrum with Spectrum viscosity Transient peak compensation step 3.16 2.93 2.95 2.9 18.38 18.36 18.37 18.2 1.83 1.75 1.76 2.1 Table 4.1: The theoretical and estimated wave speed for three different conditions. The estimated wave speed are based on three method: first, the frequency that the peak of the spectrum happens was used to reconstruct the wave speed; second, the previous value was enhanced to account for the effect of viscosity; third, the wave-front due to a transient step was analyzed to recover the speed. 97 Chapter 4. Influence of Boundary Conditions Excitation Theory (m/s) Estimation (m/s) Wavefront peak Resonance Flat plate 4.1 4.0 4.2 Point source 2.4 2.3 2.2 Table 4.2: The theoretical and estimated wave speed for an experiment with a soft PVC phantom. A flat and a point-source excitation were applied to explore the effect of boundary conditions on the speed of longitudinal waves. is delineated by a strain line passing through the peaks. A least-squares line has been fitted to the peak locations for the onset of the wave in the phantom as shown in Fig. 4.7(d). The slope of this line indicates that a longitudinal wave travels at a speed of 2.3m/s in this case. Also, from the natural vibration of the phantom, a period of 35.8ms can be infered which corresponds to a resonance frequency of 27.9Hz. Using the same calculation as for the previous case, the speed can be estimated at c = 2.2m/s. The results of these two experiments are summarized in Table 4.2. 4.6 Discussion The speed of propagation of longitudinal waves is highly dependent on the boundary conditions as well as the material properties. In nearly incompressible materials, when the boundary conditions force the wave to propagate only through a change of volume, a high speed (cp ) p-wave occurs. The p-wave, traveling as fast as ultrasound in the medium, cannot be tracked by means of ultrasound, hence compressional transient elastography has received very limited attention in the literature. However, it is shown here that with proper boundary conditions and excitation, it is feasible to induce a longitudinally propagating wave in the medium that has a much lower speed than the p-wave. Specifically, it is shown that if a block of material can freely preserve its volume, having a large uniform exciter on the surface can approximate the semi-infinite thin rod situation in which an elastic wave travels at only a few meters per second (c0 ). On the other hand, if the large exciter is replaced by a point-source indenter, the speed can be nearly reduced to the shear wave speed in the medium (c0 ). It was previously shown in [21], by using a Green’s function approach, that besides 98 Chapter 4. Influence of Boundary Conditions the s- and p-components of the wave, a coupling term is also generated in a semi-infinite medium which travels longitudinally in the near field at the speed of the shear wave. However, in this paper we elaborated that this speed depends on the dimensions of the source of excitation relative to the medium dimensions. One can expect that a mid-size exciter can produce other velocities between c0 and c0 while tuning the boundary conditions and the geometry can yield velocities between c0 and cp . In such complex cases, the velocity may not necessarily be constant throughout the medium. In a block of material, when the sides are confined, the volume change would be inevitable and therefore, it is only the p-wave that conveys the information. Due to near incompressibility of the material, the p-wave generated in the medium will be much faster than the longitudinal elastic wave when the sides are free. Finite element analysis was performed to validate the theoretical findings and explore the effect of the boundary conditions. To achieve stable finite element solutions, the Poisson’s ratio is chosen smaller than the actual value. Although a closer value to 0.5 is more realistic for incompressible soft tissues, it causes locking of the elements and instability in solving equation (4.45). Locking or hourgalssing of the elements make them behave in a stiffer fashion. The problem of locking can be reduced by numerical integration instead of closed-form integration [6]; however, this may introduce inaccuracy in the final solution. A remedy to this trade-off is to assume near incompressibility, thus ν = 0.495. As a result, locking would not happen, however the speed of the p-wave would be smaller in the simulations compared to the actual speed of ultrasound in incompressible materials where Poison’s ratio is closer to 0.5. A closer value of Poisson’s ratio to 0.5 results in higher values of λ and cp . For example, having E = 10kP a and ν = 0.5 − 7 × 10−7 gives cp = 1540m/s. The peak of the transfer functions has been used to calculate the wave speed. At the resonance frequency of a purely elastic phantom, all the points vibrate at the same phase. The first mode corresponds to the frequency that a full wavelength is trapped in the finite medium, therefore the wavelength of the wave would be twice as large as the axial dimension of the phantom. At such a frequency, the middle of the phantom experiences the highest amplitude of vibration which corresponds to a peak in the corresponding transfer function magnitude as shown in Fig. 4.4(a). The phase also drops at the natural frequency as seen in Fig. 4.4(b). The significant phase decrease in the simulations is a result of using a linear model, while real tissue often acts in nonlinear ways. As an example, human tissue and gelatin-based phantoms are known to have a power-law viscous behavior, due to which 99 Chapter 4. Influence of Boundary Conditions viscosity and relaxation-time drop as the frequency increases [8,14,15]. This nonlinearity causes the phase to have a smaller change at high frequencies. However, a linear FEM model is unable to account for such effects. As shown in Fig. 4.5, having a lower viscosity gives rise to various resonances in the system at high frequencies. Therefore, generally one may see several peaks in the transfer function, given a wide-band excitation. It is known that when a uniform quasi-static compression is applied to a block of material from one side while the other side is bounded, the axial displacements attenuate almost linearly. However, in the case of a point-source exciter, the intensity falls off similarly to an ultrasound beam as the square of the distance from the source [5]. Therefore, the quasistatic displacement amplitude for the case of a point-source is expected to be less than that of a large exciter. This explains the smaller transfer function magnitude for a point-source compared to the other two cases in Fig. 4.4(a) and the non-uniform spacing of the nodal displacement curves in Fig. 4.6(e) while Figs. 4.6(a) and 4.6(c) exhibit uniform variation between the displacements at the tail of the curves. It has been shown that a longitudinal wave may travel at a slow speed between c0 and c0 in finite medium, depending on the boundary conditions. As seen in Fig. 4.6(b) and 4.6(f), for a finite block of material a low speed elastic wave is not the only wave present in the medium. Upon exerting an axial compression to a finite block, a p-wave is first generated in the medium which travels at the high speed of cp . At this instant, the wave has not reached the boundaries and therefore, the boundary conditions and geometry are not known. Once such information is collected by the p-wave, depending on the boundary conditions, a slower wave starts to propagate that is influenced by the material properties as well as the boundary conditions. However, if the size of the block is infinite compared to the exciter such as in ultrasound or acoustic waves - bulging and rotation of infinitesimal blocks would be infeasible. In this case, the irrotational wave propagates solely through volumetric change of the medium at a speed close to cp . As noted before, when dealing with living tissue or tissue-mimicking materials where Poisson’s ratio is closer to 0.5, the speed of the p-wave is cp ≈ 1540m/s, while the speed of the slower longitudinal wave is not very sensitive to the Poisson’s ratio. However, in a bounded finite medium, the p-wave is the only type of longitudinal wave that can exist as seen in Fig. 4.6(c). Note that in this case, the resonance of the p-wave lasts for a longer time compared to the p-wave resonance in Figs. 4.6(a) and 4.6(e), since the boundary conditions assert the presence of this type of wave in the medium. 100 Chapter 4. Influence of Boundary Conditions The analytical and simulation findings were also validated experimentally on an incompressible viscoelastic phantom with known elastic properties in Fig. 4.7. Only the onset of the wave front in the medium has been used to estimate the wave speed. The reason is that the excitation was applied with a voice-coil system which could not produce a perfect step with an infinitely small rise-time due to its open-loop control. As a result, by the time the motion at the exciter has reached its peak, a wave had already started to propagate in the medium. This wave may reflect from the other side of the block and interfere with the forward-traveling wave front somewhere in the middle of the phantom depending on the speed of the wave. Due to this interference, the phase velocity changes beyond 20mm deep in the phantom in Fig. 4.7(a) and beyond 15mm in Fig. 4.7(c). How far the peak of the wavefront travels before this interference occurs depends on the rise-time of the excitation, wave speed and the length of the medium. 4.7 Conclusion It is feasible to use longitudinal waves in order to identify the mechanical properties of soft tissues. Depending on the boundary conditions, a low speed longitudinal wave can be generated in the medium which can be tracked by means of standard pulse-echo ultrasound. The speed of this wave is determined by the material properties, exciter size and other boundary conditions. The longitudinal wave speed for a viscoelastic homogeneous material can be found from the resonance frequency of the phantom, peak of the transfer function between a pair of locations, or phase speed of the transient wavefront. The analytical findings were validated numerically using 3D viscoelastic FEM models and experimentally on a tissue-mimicking phantom. 101 Chapter 4. Influence of Boundary Conditions Unbounded Elastic Excitation Unbounded Elastic Excitation Displacement Logarithm Displacement (mm) 1 0.8 0.6 0.4 0.2 0 0 20 40 60 80 Excitation Second node Center node Last node Bottom node 0 100 5 10 15 20 25 20 25 20 25 Time (ms) Time (ms) (a) (b) Bounded Medium Bounded Medium 1 Displacement (mm) Displacement (mm) 1 0.8 0.6 0.4 0.2 0.8 0.6 0.4 0.2 0 0 0 20 40 60 80 100 0 5 Time (ms) 10 15 Time (ms) (c) (d) Point−source Excitation Point−source Excitation Displacement Logaorithm Displacement (mm) 1 0.8 0.6 0.4 0.2 0 0 20 40 60 Time (ms) (e) 80 100 0 5 10 15 Time (ms) (f) Figure 4.6: Displacements at the middle line for a step excitation. The left column has a larger time-span to show the low-speed wave, where each line represents the displacement of one node. The high-speed wave can be seen in the short time-span of the right column. In (a) and (b) the excitation was applied to the top surface when the four sides were free. (c) and (d) show the response for a single-node exciter at the middle node. In (e) and (f) the medium was bounded and could not bulge. Note that the displacement axes in (b) and (d) are in a logarithmic scale. 102 Chapter 4. Influence of Boundary Conditions 20 19.0 ms Wavefront Peak Poistion (mm) Initial Position (mm) 50 40 30 20 10 0 0 50 100 Time (ms) 150 |slope| = 4.0 m/s 25 30 35 Wave travel from the plate Fitted linear model 40 13 200 14 15 16 Time (ms) (a) Wavefront Peak Position (mm) Initial Position (mm) 40 30 20 10 50 100 Time (ms) (c) 18 (b) 35.8 ms 0 0 17 150 200 Wave travel from the point−source Fitted linear model 30 35 |slope| = 2.3 m/s 40 11 12 13 Time (ms) 14 15 (d) Figure 4.7: Wave propagation in a soft PVC phantom with E = 17kP a and an axial length of 4cm under two different boundary conditions when a step excitation is applied. (a) and (c) show the displacements of tissue at different depths for a flat and a point-source excitation respectively. The location of the wavefront peaks are shown versus time for these two boundary conditions in (b) and (d), where a line has been fitted to the data to estimate the speed. The resonance frequency of the phantom is also a function of the elastic modulus, geometry and boundary conditions. 103 References [1] Baghani, A., Eskandari, H., Salcudean, S. E. and Rohling, R.: 2009, Measurement of Viscoelastic Properties of Tissue Mimicking Material Using Longitudinal Wave Excitation, Submitted to the IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control. [2] Bercoff, J., Tanter, M. and Fink, M.: 2004, Supersonic shear imaging: a new technique for soft tissue elasticity mapping, IEEE Trans Ultrasonics, Ferroelectrics and Frequency Control 51(4), 396–409. [3] Bro-Nielsen, M.: 1998, Finite element modeling in surgery simulation, Proceedings of the IEEE 86(3), 490–503. [4] Catheline, S., Gennisson, J., Delon, G., Fink, M., Sinkus, R., Abouelkaram, S. and Culioli, J.: 2004, Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: An inverse problem approach, Journal of the Acoustical Society of America 116(6), 3734–3741. [5] Cobbold, R. S. C.: 2007, Foundations of Biomedical Ultrasound, Oxford University Press, NY. [6] Cook, R. D., Malkus, D. S. and Palesha, M. E.: 1989, Concepts and Applications of Finite Element Analysis, 3rd edn, John Wiley & Sons, Inc., New York. [7] Eskandari, H., Salcudean, S. E. and Rohling, R.: 2007, Tissue Strain Imaging Using a Wavelet Transform-Based Peak Search Algorithm, IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control 54(6), 1118-1130. [8] Eskandari, H., Salcudean, S. E. and Rohling, R.: 2008a, Viscoelastic parameter estimation based on spectral analysis, IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control 55(7), 1611–1625. 104 Chapter 4. Influence of Boundary Conditions [9] Eskandari, H., Salcudean, S. E., Rohling, R. and Ohayon, J.: 2008b, Viscoelastic characterization of soft tissue from dynamic finite element models, Physics in Medicine and Biology 53(22), 6569-6590. [10] Fatemi, M. and Greenleaf, J. F.: 1999, Vibro-acoustography: An imaging modality based on ultrasound-stimulated acoustic emission, Proceedings of the National Academy of Science, U.S.A. 96(12), 6603– 6608. [11] Fu, D., Levinson, S. F., Gracewski, S. M. and Parker, K. J.: 2000, Non-invasive quantitative reconstruction of tissue elasticity using an iterative forward approach, Physics in Medicine and Biology 45(6), 1495– 509. [12] Insana, M. F., Pellot-Barakat, C., Sridhar, M. and Lindfors, K. K.: 2004, Viscoelastic imaging of breast tumor microenvironment with ultrasound, Journal of Mammary Gland Biology and Neoplasia 9(4), 393– 404. [13] Insana, M., Liu, J., Sridhar, M. and Pellot-Barakat, C.: 2005, Ultrasonic mechanical relaxation imaging and the material science of breast cancer, Ultrasonics Symposium, 2005 IEEE, Vol. 2, pp. 739–742. [14] Joly-Duhamel, C., Hellio, D. and Djabourov, M.: 2002, All gelatin networks: 1. Biodiversity and physical chemistry, Langmuir 18(19), 7208– 7217. [15] Kiss, M., Varghese, T. and Hall, T.: 2004, Viscoelastic characterization of in vitro canine tissue, Physics in Medicine and Biology 49(18), 4207– 18. [16] Kolsky, H.: 1963, Stress Waves in Solids, Dover Publications Inc., NY. [17] Ljung, L.: 1999, System identification : theory for the user, Prentice Hall PTR. [18] Malvern, L. E.: 1969, Introduction to the mechanics of a continuous medium, 1st edn, Prentice-Hall Inc., New Jersey. [19] McLaughlin, J. and Renzi, D.: 2006, Shear wave speed recovery in transient elastography and supersonic imaging using propagating fronts, Inverse Problems 22(2), 681–706. 105 Chapter 4. Influence of Boundary Conditions [20] Nightingale, K., McAleavey, S. and Trahey, G.: 2003, Shear-wave generation using acoustic radiation force: in vivo and ex vivo results, Ultrasound in Medicine and Biology 29(12), 1715–1723. [21] Sandrin, L., Cassereau, D. and Fink, M.: 2004, The role of the coupling term in transient elastography, J Acoust Soc Am 115(1), 73–83. [22] Sandrin, L., Tanter, M., Catheline, S. and Fink, M.: 2001, Timeresolved 2d pulsed elastography: experiments on tissue-equivalent phantoms and breast in vivo, Proceedings of SPIE, Vol. 4325, pp. 120– 126. [23] Sandrin, L., Tanter, M., Gennisson, J.-L., Catheline, S. and Fink, M.: 2002, Shear elasticity probe for soft tissues with 1-d transient elastography, IEEE Trans Ultrasonics Ferroelectrics and Frequency Control 49(4), 436–446. [24] Sarvazyan, A. P., Rudenko, O. V., Swanson, S. D., Fowlkes, J. B. and Emelianov, S. Y.: 1998, Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics, Ultrasound in Medicine and Biology 24(9), 1419–1435. [25] Sinkus, R., Tanter, M., Xydeas, T., Catheline, S., Bercoff, J. and Fink, M.: 2005, Viscoelastic shear properties of in vivo breast lesions measured by MR elastography, Magnetic Resonance Imaging 23(2), 159–65. [26] Thompson, J. H. C.: 1933, On the theory of visco-elasticity: A thermodynamical treatment of visco-elasticity, and some problems of the vibrations of visco-elastic solids, Philosophical Transactions A 231, 339–407. [27] Turgay, E., Salcudean, S. E. and Rohling, R.: 2006, Identifying the mechanical properties of tissue by ultrasound strain imaging, Ultrasound in Medicine and Biology 32(2), 221–235. [28] Zahiri-Azar, R. and Salcudean, S.: 2006, Motion estimation in ultrasound images using time domain cross correlation with prior estimates, IEEE Transactions on Biomedical Engineering 53(10), 1990–2000. [29] Zhang, M., Castaneda, B., Wu, Z., Nigwekar, P., Joseph, J. V., Rubens, D. J. and Parker, K. J.: 2007, Congruence of imaging estimators and mechanical measurements of viscoelastic properties of soft tissues, Ultrasound in Medicine and Biology 33(10), 1617–1631. 106 Chapter 5 A High Frame Rate Ultrasound System for the Study of Tissue Motions∗ 5.1 Introduction Medical ultrasound imaging systems are used to image both the echogenicity of the tissue, as in B-mode imaging, and the tissue motion, as in Doppler imaging [1, 2]. Over the past two decades considerable advances have been made in the field of tissue motion imaging. State-of-the-art speckle tracking software is capable of measuring sub-micron displacements inside tissue [3– 11]. A fundamental limitation of ultrasound based motion tracking is the inherent low frame rate of the ultrasound systems. Due to the finite speed of the ultrasound waves in human tissue, acquiring a full 2D B-mode image requires tens of milliseconds. A low frame rate limits the scope of the motions that can be properly detected by the system; for a typical frame rate of 100Hz, the Nyquist sampling criterion implies that for a proper estimation, the tissue motions should be band limited to 50Hz. Not all of the natural body motions meet this criterion. In cardiovascular imaging, higher frame rates are required [2]. Motions induced inside the tissue to study its mechanical characteristics, as in dynamic elastography, surpass this limit in many cases [12–16]. Different techniques have been devised over the years to increase the frame rate of ultrasound systems for motion estimation. In conventional pulsed Doppler imaging, a single scan-line is acquired repeatedly at a high ∗ A version of this chapter has been submitted for publication: A. Baghani, A. Brant, S.E. Salcudean, R. Rohling, “A High Frame Rate Ultrasound System for the Study of Tissue Motions.” 107 Chapter 5. High Frame Rate Ultrasound frame rate, the so called pulse repetition frequency (PRF) in that context [1, 2]. This frequency also determines the range of the speeds which the system can measure. Higher speeds demand a higher pulse repetition frequency. In conventional color (and power) Doppler imaging, the window of interest is subdivided into sectors [1, 2]. The width of the sectors is determined by the desired PRF, which is set by the user. Each sector is imaged repeatedly a number of times at the PRF, and velocity estimations are carried out. Then the system moves on to acquiring the next sector. When all the sectors have been acquired, the velocity data for the whole window of interest is updated on the screen. Sonoelastography directly uses this data for elasticity imaging [12]. The idea of dividing the region of interest into a number of sectors, each of which can be acquired separately at a higher frame rate, has been used in a number of high frame rate imaging systems. Synchronization of the excitation with the sector acquisition is used in [17, 18] to achieve a frame rate of 3kHz. Synchronization of the sector acquisition with the heart beat through echocardiogram (ECG)-gating is used in [19,20] to achieve an effective frame rate of 481Hz at an imaging depth of 11cm and 64 line density for cardiac imaging. In [21], the authors use morphing to produce additional frames between successive acquired frames. The velocity data obtained from the speckle tracking is used in the morphing algorithm. Although this method helps the clinicians by increasing the frame rate of the video seen on the screen, it does not overcome the inherent Nyquist limit for measuring high frequency components of the tissue motion. Another approach to high frame rate ultrasound imaging is to use nonconventional pulsing techniques. Ultrafast ultrasound is one such technique in which transmit focusing is not applied [22]; all the probe crystals are fired at once to create a linear wave front. Dynamic receive focusing is applied to all the received signals by post processing to form the equivalent RF-lines for each crystal. The lack of transmit focusing slightly lowers the quality of B-mode images. However very good results have been obtained in detecting high frequency motions. The system has been implemented in AixplorerTM , a next generation ultrasound machine marketed by SuperSonic Imagine Co., Aix-en-Provence, France. Dedicated hardware for parallel acquisition and processing of data from an aperture of crystals is required in these systems [23]. A number of coded pulsing techniques have also been proposed which can increase the frame rate [24, 25]. Frequency modulation was used to 108 Chapter 5. High Frame Rate Ultrasound increase the frame rate by a factor of 5 and Hadamard spatial encoding in combination with frequency modulation increased the frame rate by a factor of 25 [24]. A slight degradation in the B-mode image quality is reported with these methods. However, the performance of the system in motion tracking has not been studied. In [25] binary codes are used on the transmitted signals to produce distinct impulse responses from different directions in the region of interest. The received signals are then processed in parallel to decode the data and form the image. Since multiple receive directions are compounded, blocking artifacts occur. The performance of the system in motion tracking has not been studied. In this article we report a high frame rate ultrasound system which uses the concept of sector subdivision and can readily be implemented on conventional ultrasound systems. A detailed analysis is carried out and conditions set forth on the tissue motions that can accurately be estimated using this system. The main novelty of the system is in its phase correction scheme which makes it possible to make virtually simultaneous motion measurements at all points throughout the tissue being imaged. Another advantage of the system is its real-time capability. This enables the user to overlay motion estimation images on B-mode images while the probe is scanning. The system has been implemented on a Sonix RP machine (Ultrasonix Medical Corporation, Richmond, BC, Canada). Examples of the applications of the system are given. The rest of this article is organized as follows. In Section 5.2 the high frame rate acquisition technique is reviewed and the novel phase correction scheme presented. Section 5.3 presents a set of conditions which should be satisfied by the tissue motion. These conditions ensure that the system estimates the tissue motion correctly. Section 5.4 presents some applications of the system with experimental results. We conclude the article in Section 5.5. 5.2 5.2.1 Sector-Based High Frame Rate Acquisition Technique Sequencing In this section we discuss the sector-based sequencing scheme on which the high frame rate data acquisition is based. In this discussion, instead of the term “firing an element” we use the term “firing a scan-line” by which we 109 Chapter 5. High Frame Rate Ultrasound mean that the usual beam forming technique for focusing the beam at a certain depth is performed and an aperture of elements are fired. In conventional B-mode imaging, the scan-lines are fired one by one in a sequence starting from the first line and finishing by the last line. This constitutes a single frame of the data. Figure 5.1 shows the B-mode acquisition time-line for two frames on a twelve element probe. Each slanted line represents the points at different depths on one of the scan-lines. The number of the scan-line is written on the slanted line. Figure 5.1: Acquisition time-line for two frames of B-mode data on a twelve element probe. Consider a typical point at depth d on the typical scan-line No. 9. This point is imaged at time t1 in the first frame and at time t2 in the second frame. If we denote the time required to acquire a single line by T¯, the time interval between successive acquisitions of the point is, t2 − t1 = Ts = 12T¯. (5.1) In motion estimation algorithms, the successive acquisitions of the RFdata from each point are used to estimate the motion at that point. Hence the displacements are sampled at the frequency, fs = 1 1 1 = Ts 12 T¯ (5.2) This frequency is the same for all the points on all the twelve scan-lines of 110 Chapter 5. High Frame Rate Ultrasound the probe in this sequencing scheme, and is equal to the frame rate of the B-mode. For a typical probe of 128 elements and a typical imaging depth of 6cm this frequency is equal to 100Hz. Figure 5.2 shows a typical sector-based sequencing with the twelve element probe which would increase the Nyquist limit by a factor of 3. In this sequence the first four firings are the same as the B-mode imaging, i.e. firing scan-lines one to four. However instead of moving on to scan-line five in the fifth firing, the sequencer goes back to the first scan-line and repeats the first four scan-lines a couple of times (just once in this case for illustration) before moving on to the fifth scan-line. The twelve scan-lines of the probe are thus divided into three sectors (each, four scan-lines wide), and each sector is scanned twice (a multitude of times in general) before the next sector is acquired. Figure 5.2: Acquisition time-line for two frames of high frame rate acquisition data on a twelve element probe As can be seen from Fig. 5.2, the time interval between successive scannings of the typical point at depth d on scan-line No. 9 has been reduced. In other words the sampling frequency has been increased during the interval the point is being observed (two frames in this example for illustration). 111 Chapter 5. High Frame Rate Ultrasound The sampling frequency in this case is, fs = 1 11 = ¯ Ts 4T (5.3) which is 3 times higher than that achieved in (5.2). In general, the sector size can be chosen arbitrarily. If M scan-lines are grouped together in a sector and the sector is scanned K times in a sequence, the sampling frequency would be, fs (M ) = 1 1 M T¯ (5.4) From this formula, the highest frame rate is achieved with M = 1. In this case each line is scanned K times, before scanning the next line. Note that the total time needed to acquire K complete frames does not depend on M , the number of scan-lines grouped together in a sector, but on N , the total number of crystals in the probe, Ttotal = K(M T¯) N = KN T¯ M (5.5) Note that for the example of Fig. 5.2, the number of lines in each sector M is equal to 4, the number of repetitions (frames) K is equal to 2, and the total number of elements in the probe N is equal to 12. The number of sectors needed to cover the entire probe is, N 12 = = 3. M 4 (5.6) As noted in the introduction, the presented sector-based sequencing scheme is used in conventional color and power Doppler imaging. In these imaging modalities a window (region of interest) is used which reduces the total number of lines scanned, N . The number of lines in each sector, M , is determined based on the PRF which is itself determined by the desired speed range. This frequency is the same as fs (M ) presented in (5.4). The presented sequencing scheme does not change the overall acquisition speed in terms of the total time needed for acquiring a number of frames. However it increases the sampling frequency by changing the order in which the lines are acquired. This reordering also changes the delay pattern between the acquisition of the scan lines. In the next subsection we present our novel phase correction algorithm for compensating these delays. 112 Chapter 5. High Frame Rate Ultrasound 5.2.2 Time Delay Compensation The theory of discrete signals and systems deals with discretizations of continuous-time signals and systems obtained through sampling. By far, the most common type of sampling is the one in which all signals and systems are sampled simultaneously and at a fixed rate. Although cases where the sampling rate changes or the signals are not sampled simultaneously are less common, the theory still provides the necessary tools for dealing with them. Ultrasound based motion tracking is one. Tissue which is undergoing dynamic deformations can be considered to be a continuous-time dynamic system. The signals of interest are the positions or the velocities of the different points. It is of interest to sample all these signals simultaneously. An ideal time-line for such an acquisition is shown in Fig. 5.3. A comparison with Fig. 5.1 shows that, • the scan-lines are no longer slanted, which means that all the points at all depths are acquired instantly at the tick of the sampling clock, and • the scan-lines are not separated in time, which means that all the scan-lines are acquired simultaneously. In reality, because of the finite propagation speed of the ultrasound inside tissue, this is impossible to achieve. The time delay from the tick of the sampling clock to the acquisition time of a point depends on the sequencing method used and the position of the point in the sequence. Figure 5.3: A hypothetical acquisition time-line for delay-free simultaneous sampling on a twelve element probe 113 Chapter 5. High Frame Rate Ultrasound For periodic motions, the time delays can be compensated in the frequency domain. Assume that, • the probe does not move with respect to the tissue, • the speed of ultrasound is constant throughout the tissue and equal to c, • and that all the points to be scanned are vibrating at a single frequency of fe . For motions containing a multitude of frequencies, the process to be presented can be repeated for each frequency separately. The acquisition process is triggered by the sampling clock at time zero. We denote the measured displacement phasor (at frequency fe ) of the points on line m in sector n by, Umeas (d; m, n) m = 1, . . . , M n = 1, . . . , N/M (5.7) where d, as before, denotes the depth on the line. We start with the time delays in a single scan-line. Since the ultrasound has a finite traveling speed c, the points on a single line are not all scanned at the same instant. A point at depth d inside the tissue would be scanned t = 2d/c later than the first point on the line. Since the point is vibrating at frequency fe , by this time, its phase has advanced by, φ = 2πfe t = 2πfe 2d . c (5.8) The phasor therefore needs to be compensated as, Ucomp1 (d; m, n) = Umeas (d; m, n) exp −j2πfe 2d c (5.9) The next time delay to deal with is the time delay from the acquisition of the first scan-line to the scan-line m in a sector. From Fig. 5.2, the delay between the two lines is (m − 1)T¯. Therefore, the phasor needs to be compensated as, Ucomp2 (d; m, n) = Ucomp1 (d; m, n) exp −j2πfe (m − 1)T¯ (5.10) The last time delay to deal with is the time delay from the acquisition of the first scan-line in sector one to the first scan-line in sector n. For M lines in a sector, each scanned K times, this time delay is (n − 1)M K T¯. 114 Chapter 5. High Frame Rate Ultrasound Therefore the phasor needs to be compensated as, Ucomp3 (d; m, n) = Ucomp2 (d; m, n) exp −j2πfe (n − 1)M K T¯ (5.11) The phasor Ucomp3 (d; m, n) contains the amplitude and the corrected phase information for all the points, which by assumption are all vibrating at frequency fe . This phasor could be used in the frequency domain for inverse problem calculations, or converted into the time domain, u(t; m, n) = |Ucomp3 (d; m, n)| cos (2πfe t + ∠Ucomp3 (d; m, n)) (5.12) To have actual numbers for these delays we consider a simple example. Assume that the probe has 128 elements, N = 128, the sectors are 16 elements wide, M = 16, and 50 frames are acquired, K = 50. Assume that the ultrasound speed is 1540m/s and that the depth of imaging is 77mm. In this case the time required to acquire a single line is, 2 × 77mm T¯ = = 100µs 1540m/s (5.13) Assume that the tissue is vibrating at a frequency of 100Hz, fe = 100Hz. For each line the required phase compensation is in the range, φ1 ∈ [0 , 2π × 100Hz × 100µs] = [0 , 0.0628] . (5.14) The required phase compensations for the sixteen lines in each sector are, φ2 ∈ {0.0628 × (m − 1)|m = 1, . . . , 16} ∈ {0, 0.0628, 0.1257, . . . , 0.9425} . (5.15) The required phase compensations for different sectors are, φ3 ∈ {0.0628 × 16 × 50 × (n − 1)|n = 1, . . . , 8} ∈ {0, 50.2655, 100.531, . . . , 351.8584} . (5.16) From (5.5), the total time required to acquire the entire data set is, Ttotal = 50 × 128 × 100µs = 0.64s . (5.17) 115 Chapter 5. High Frame Rate Ultrasound 5.2.3 Implementation Issues The scenario presented in the previous subsection assumes that the ultrasound imaging system could be programmed for a sequence of any length. In practice, there is a limitation on the number of scan-lines that can be programmed into a sequencer. The sequencer then loops through that sequence repeatedly. For the Ultrasonix RP system which is used for this research, a typical limit of 500 scan-lines could be considered (the exact number depends on different settings such as depth and sampling frequency). For a 128 element probe, if K = 100 samples (frames) are required, the total number of lines to be programmed, if the entire sequence is to be programmed in one batch, is equal to 12800 scan-lines. This number is far beyond the 500 limit. The solution is to program each sector individually and let the sequencer loop through that sector the desired number of frames, K, and then reprogram the sequencer for the next sector. For instance if the sector size is M = 8, and K = 100 samples are required, this method requires only 8 lines to be programmed and the sequencer repeats the sequence 100 times. The sequencer is then stopped and the next sector programmed. A flowchart of the process is shown in Fig. 5.4. The drawback of this solution is that there will be a (variable) time delay between the acquisitions of the different sectors. This is due to the time required for the software to download the new sector to the sequencer hardware. A typical time-line for this acquisition scheme is presented in Fig. 5.5. The variable delay invalidates the third compensation step carried out in the previous subsection. The time delay for inter-sector compensation is not (n − 1)M K T¯ anymore, but an unknown value. Different solutions are available to this problem. One approach is to use time stamps. If the time that the data collection is started for each sector is saved, say in Ti , the phasors can be compensated as, Ucomp3 (d; m, n) = Ucomp2 (d; m, n) exp (−j2πfe (Tn − T1 )) (5.18) Another approach is to choose the sectors so that they overlap. For instance if the first sector consists of scan-lines 1 to 4, the second sector can consist of scan-lines 4 to 7. By selecting the sectors in this way, the scanline No. 4 is common to both sectors. Now there are two measurements available for this line, which after the first two steps of compensation become 116 Chapter 5. High Frame Rate Ultrasound Figure 5.4: Flow chart of sequence programming for high frame rate image acquisition 117 Chapter 5. High Frame Rate Ultrasound Figure 5.5: Acquisition time-line involving reprogramming delay on a twelve element probe Ucomp2 (d; 4, 1) and Ucomp2 (d; 1, 2). Since we have assumed that the probe is not moving with respect to the tissue during the acquisition process, these two phasor profiles should be exactly the same for all depths d, except for a constant phase delay. The phase delay can be found by a simple linear least squares: min Ucomp2 (d; 1, 2) − (αUcomp2 (d; 4, 1) + β) 2 . (5.19) α,β∈C The third compensation step becomes Ucomp3 (d; m, 2) = αUcomp2 (d; m, 2) . (5.20) This procedure is repeated for all the inter-sector delays. Another approach is to make predictions. Assume that the first sector consists of scan-lines 1 to 4 and the second sector consists of scan-lines 5 to 8. Based on scan-lines 1 to 4 a prediction can be made about the phasor of the points on scan-line No. 5. For instance a linear extrapolation in time from scan lines 3 and 4 yields, Upredicted (d; 1, 2) = 2Ucomp2 (d; 4, 1) −Ucomp2 (d; 3, 1) . (5.21) 118 Chapter 5. High Frame Rate Ultrasound This predicted phasor profile and the measured phasor profile after compensation, Ucomp2 (d; 1, 2), should be almost the same for all depths d, except for a constant phase delay. The phase delay can be found and compensated for by a similar procedure to (5.19) and (5.20). This is the solution we have implemented in our system. 5.3 Requirements on the Tissue Motion As was shown in the previous section, if all the points are vibrating at a frequency of fe in the steady state, it is possible to compensate for the acquisition phase lags in the phasors. Now we address the question of whether it is possible to obtain the phasors of displacements for every frequency fe . To further elaborate on the question, consider a single point inside the tissue which is vibrating at a frequency of fe . The displacement of this point is sampled at a frequency of fs (M ) given by (5.4). Following the Nyquist criterion, to be able to recover the correct amplitude and phase of the displacement, it is necessary that, 1 fe < fs (M ). 2 (5.22) If it were possible to observe the displacement for an infinite amount of time and obtain infinitely many samples, (5.22) would be the only restriction on the tissue motion. However a review of Subsection 5.2.1 reveals that for each point, only K samples are taken, and the observation time span for each point is, K Tobs = . (5.23) fs (M ) This has very important consequences on the applicable values of fe . With reference to Fig. 5.5 each sector is being observed during a different interval in time. Conceptually speaking, if the collected data from the different sectors are to be brought together and reconciled, the motions should be periodic with the observation interval being an integer multiple of the full period. In Fourier theory, a signal with a limited time duration Tobs is analyzed with a Fourier series. The implicit assumptions in the conversion of such a signal to the frequency domain is that the signal is part of a periodic signal, and that the observation interval Tobs is a full period of the periodic signal. Figure 5.6 shows the extension of a limited duration signal to a periodic signal. 119 Chapter 5. High Frame Rate Ultrasound Figure 5.6: The extension of a limited duration signal to a periodic signal Figure 5.7: Signal processing pipeline of the high frame rate ultrasound system 120 Chapter 5. High Frame Rate Ultrasound In the frequency domain, the power spectrum of a periodic signal y(t) of period Tobs contains energy only at integer multiples of the fundamental frequency, y(t) = al exp(j2πlffund ) , (5.24) l∈Z ffund = 1 Tobs . (5.25) Therefore the motion should only have components at frequencies, fe = lffund = lfs (M ) K l ∈ Z. (5.26) It is worth mentioning that this type of discrete spectrum by no means limits the generality of the shape of the signal over the time duration Tobs . However the Nyquist criterion (5.22) limits the range of the permissible frequencies for the tissue motion. Substitution of (5.26) into the criterion gives − K K <l< . 2 2 (5.27) In summary, the tissue motion needs to be a periodic band-limited motion. Most of the natural body motions such as the motions resulting from the heartbeat are periodic. However the period might change slightly between different pulses. A method for compensating these variations can be found in [20]. If the imaging system is used together with an excitation technique for elastography applications, the total number of frequencies which could be used to excite the tissue is equal to K, the total number of times each point (or scan-line) is scanned in the sequence. The excitation frequency is always a multiple of the fundamental frequency, the inverse of the observation period. However there is no need for the excitation to be of a single frequency (harmonic) nature. As a matter of fact any combination of the K permissible frequencies could be used. The best combination may depend on the particular application. For our experiments, custom software was developed that enables the user to choose the line duration T¯, the sector width M , and the number of samples K. It computes the permissible excitation frequencies and presents them as a list from which the user can choose the frequencies for which the phasors are computed. In between the high frame rate acquisition sequences, B-mode scans are carried out. The software can then overlay the amplitude 121 Chapter 5. High Frame Rate Ultrasound or the phase of the computed phasors on the B-mode images in real-time. The latency depends on the total acquisition time given in (5.5). For a typical probe of N = 128 crystals and a typical imaging depth of d = 6cm if K = 100 samples are taken the latency is 1 second. The signal processing pipeline is shown in Fig. 5.7. 5.4 Applications and Experimental Results In this section applications of the high frame rate system in dynamic elastography are presented. In ultrasound elastography, the tissue is subjected to a constant or time-varying force and the response of the tissue to this force is monitored by an ultrasound system. The response consists of the tissue motions at different locations as a function of time. The developed high frame rate system makes it possible to monitor tissue motions for excitations covering a wider range of frequencies. Moreover, since most of the motion-tracking algorithms use a similarity measure between two frames to estimate the motion, at higher frame rates (shorter intervals) the similarities are easier to detect, and hence these algorithms generally perform better. To test the high frame rate system, an experimental setup is developed (Fig. 5.8). The setup consists of a small table with an opening through which the ultrasound probe is fitted from below. The test sample is placed on top of the table over the ultrasound probe. The exciter consists of a speaker and a sound tube. The sound tube acts as a wave guide between the speaker and is made of plexi-glass. The inner and outer diameters are 57mm and 64mm, respectively. The tube is clamped vertically to a column and its height can be adjusted so that the open end of it touches the top surface of the test sample. The pressure variations caused by the sound waves inside the tube create a dynamic force excitation on the surface of the sample. Note the choice of the axes in Fig. 5.8; the axial, lateral and elevational directions in the ultrasound image are along the x, y and z axes respectively. In this coordinate system, the depth d is the axial coordinate x, and the sector index variables m and n determine the lateral coordinate y: U (d; m, n) → U (x, y) (5.28) The first test is carried out on a cylindrical phantom. The phantom was molded out of Super Soft PlasticTM (M-F Manufacturing Co., Inc. Fort Worth, Tx, USA), a PolyVinyl Chloride (PVC)-based plastic. SigmaCell Cellulose, Type 50 (Sigma-Aldrich Inc., St Louis, MO, USA) was added as 122 Chapter 5. High Frame Rate Ultrasound Figure 5.8: Schematics of (a) the side view and (b) the front view of the experimental setup developed for testing the system. Note the choice of the coordinate axes. 123 Chapter 5. High Frame Rate Ultrasound Table 5.1: Parameters of the high frame rate system parameter value probe type L9–4 linear array center frequency 6.5MHz sampling frequency 20MHz depth of focus 30mm imaging depth 56mm total number of crystals (N ) 128 sector width (M ) 8 line duration (T¯) 100µs number of samples or frames (K) 50 the ultrasound scatterer to develop the speckle pattern required for motion tracking. The diameter of the cylindrical phantom is 90mm and its height is 65mm. The parameters of the high frame rate system are summarized in Table 5.1. In this case the sampling frequency is 1.25kHz. Based on the equations (5.26) and (5.27) the permissible excitation frequencies are 0, 25, 50, 75, 100, . . . , and 625Hz. In the first experiment, a sinusoidal excitation at a frequency of 100Hz is used. This excitation is created using MATLAB and applied to the speaker through the output of the sound card on the PC. The speaker is an offthe-shelf PC-speaker and comes with an amplifier. The axial displacements D(x, y, k) are calculated by a correlation-based speckle tracking software [9]. Note that k is the discrete time or the frame number. The phasor of the displacements at the frequency of excitation is then computed, k=50 exp −j2π U (x, y) = k=1 k · 100Hz 1250Hz D(x, y, k) . (5.29) The phase of the phasor is then compensated to account for the delays. Fig. 5.9 illustrates the steps carried out for the phase compensation. In Fig. 5.9(a) The original phase of U (x, y) is plotted. The sectors, which are 8 scan-lines wide, can clearly be distinguished. Fig. 5.9(b) shows the effect of the phase compensation for the delays on each scan line, i.e. the phase of Ucomp1 (x, y). Since the compensated phase is small, the difference is not noticeable. Fig. 5.9(c) shows the effect of the phase compensation for the delays between the scan lines in each sector, i.e. the phase of Ucomp2 (x, y). The difference is now more evident, as the slope of the phase changes can 124 Chapter 5. High Frame Rate Ultrasound be seen clearly to have been reduced. Fig. 5.9(d) shows the effect of the phase compensation for the delays between the sectors, i.e. the phase of Ucomp3 (x, y). The phase is clearly smooth in this figure. This shows the effectiveness of the phase correction method. Figure 5.9: Different steps in the phase correction process: the phase (a) before any compensation (b) after compensation for the delays on each scanline (c) after compensation for the delays between the scan-lines in each sector (d) after compensation for the delays between the sectors. In the second experiment, a multiple frequency excitation is used on the same phantom. The excitation used is of the following form: excitation = 1 sin(2π75Hz t) + 1 sin(2π100Hz t) + 2 sin(2π150Hz t) + 3 sin(2π200Hz t) (5.30) Figure 5.10 shows the amplitude |U (x, y)| and corrected phase ∠U (x, y) images of the axial displacement phasors at the frequencies of excitation. The white lines in the phase images are contours of zero phase. As can be seen all images show smooth variations. The nodes and peaks become more frequent with the increase in frequency. Note that without the phase correction, the amplitude images would still be smooth. However, the phase 125 Chapter 5. High Frame Rate Ultrasound images would be striped as in Fig. 5.9(a). Figure 5.10: (a), (b), (c), (d) the amplitude, and (e), (f), (g), (h) corresponding phase of the displacement phasors at different frequencies of excitation. From left to right the frequencies are 75, 100, 150, and 200Hz. The phasor image at 100Hz was chosen, and its spatial derivatives calculated. The results are presented in Fig. 5.11. Note that without the phase correction it would not be possible to compute the lateral derivatives, ∂U/∂y and ∂ 2 U/∂y 2 . Because of the delay in each scan-line, the axial derivatives ∂U/∂x and ∂ 2 U/∂x2 , would also contain slight errors. The spatial derivatives of the displacement field are used extensively in elastography for solving the inverse problem. The axial strain is given by ∂U/∂x. For calculating the rotation in the xy-plane, ω ¯ z , ∂U/∂y is needed: 2¯ ωz = ∂U ∂V + , ∂y ∂x (5.31) where V (x, y) is the lateral displacement. The Laplacian of U (x, y, z) is usually approximated by: ∇2 U (x, y, z) = ≈ ≈ ∂2U ∂2U ∂2U + + ∂x2 ∂y 2 ∂z 2 2 2 ∂ U ∂ U + 2 ∂x ∂y 2 ∂2U . ∂x2 (5.32) (5.33) 126 Chapter 5. High Frame Rate Ultrasound Figure 5.11: Spatial derivatives of the axial displacement phasor U (x, y) at the excitation frequency of 100Hz. (a) |∂U/∂x| (b) |∂ 2 U/∂x2 | (c) |∂U/∂y| (d) |∂ 2 U/∂y 2 | (e) ∠∂U/∂x (f) ∠∂ 2 U/∂x2 (g) ∠∂U/∂y (h) ∠∂ 2 U/∂y 2 To demonstrate the effectiveness of the system for dynamic elastography, two inclusion phantoms were tested and their elastograms obtained. We chose direct inversion of the wave equation for obtaining the elastograms, E(x, y) = ρ (j2πfe )2 U (x, y) ∂2 ∂x2 + ∂2 ∂y 2 , (5.34) U (x, y) where fe , as before, is the frequency at which the phasor is computed (the excitation frequency). Note that the factor j2πfe in the numerator is the frequency domain equivalent of partial derivative with respect to time t. Both phantoms are rectangular with dimensions 30mm×55mm×65mm. Each phantom includes a cylindrical inclusion with a diameter of 15mm. For one of the phantoms the inclusion is softer, and for the other one it is harder than the surrounding. The blocks and the inclusions were made from Super Soft Plastic with different concentrations of Plastic Softener or Plastic Hardener. A schematic of the phantom geometry and imaging plane is presented in Fig. 5.12. The setup of Fig. 5.8 is used to excite and image the phantoms. The 127 Chapter 5. High Frame Rate Ultrasound Figure 5.12: Outline of the geometry of the inclusion phantoms used in the elastacity experiments. excitation used is of the following form: excitation = 1 sin(2π100Hz t) + 1.5 sin(2π125Hz t) − 1.5 sin(2π150Hz t) − 2 sin(2π175Hz t) + 2 sin(2π200Hz t) . (5.35) For each of the excitation frequencies, the phasor U (x, y) was calculated and was used in (5.34) to obtain an elastogram E(x, y). The obtained elastograms for different frequencies were then averaged to obtain the final elastogram. The results for the hard and soft inclusion phantoms are presented in Fig. 5.13 and 5.14 respectively. Since the echogenicity of the inclusions are close to their surrounding material, it is difficult to recognize them in the Bmode images. However, they can clearly be recognized from the elastograms in both cases. 5.5 Summary and Conclusion A high frame rate ultrasound system is introduced in this article which can readily be implemented on conventional ultrasound systems without 128 Chapter 5. High Frame Rate Ultrasound Figure 5.13: (a) B-mode and (b) elastogram images of a phantom with a circular hard inclusion. Figure 5.14: (a) B-mode and (b) elastogram images of a phantom with a circular soft inclusion. 129 Chapter 5. High Frame Rate Ultrasound the need for any hardware modifications. The system uses the concept of sector subdivision for sequencing which is widely used in Doppler imaging, to achieve high frame rates. For a sector of width one, i.e. a single scan-line, the frame rate is only limited by the time it takes for the deepest point to be scanned. A phase correction scheme was presented to compensate for the time delays involved in the acquisition process. Necessary conditions on the frequency content of the tissue motion were presented, in order for the displacements to be recovered properly. Since the system simply reorders the pulse sequence, the tissue is still imaged one line at a time. As a result the total acquisition time with this system can be a fraction of a second. The probe and tissue should be stationary during this acquisition interval to prevent motion artifacts. The system can be used together with periodic excitations to study tissue characteristics, and as such is not meant as a replacement for other high speed imaging techniques used in applications with different requirements. The phase compensation scheme presented here can readily be extended to the 3D case. For 2D phased array transducers the extension is straight forward. For mechanically actuated 3D probes with linear phased array transducers, the motion of the stepper motor inside the probe can be synchronized with the sector acquisition to produce known delays between different image planes. These known delays can then be compensated for in a similar manner. The high frame rate system can be used together with the exciters used in elastography techniques, such as mechanical exciters and acoustic radiation force (ARF) exciters, to obtain displacement data of the tissue at different frequencies of excitation. The excitation frequency can go up to a few kiloHertz and the system can still estimate the displacement data. These data can then be used in the inverse problem solving techniques to recover the tissue elasticity and viscosity. 130 References [1] R. Cobbold, Foundations of Biomedical Ultrasound. NY: Oxford University Press, 2007. [2] J. Jensen, Estimation of Blood Velocities Using Ultrasound: A Signal Processing Approach. Cambridge University Press, 1996. [3] I. Hein and W. O’Brien, “Current time-domain methods for assessing tissue motion by analysis from reflected ultrasound echoesa review,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 40, pp. 84–102, 1993. [4] E. E. Konofagou and J. Ophir, “Precision estimation and imaging of normal and shear components of the 3d strain tensor in elastography,” Phys. Med. Biol., vol. 45, pp. 1553–1563, 2000. [5] R. Righetti, J. Ophir, and P. Ktonas, “Axial resolution in elastography,” Ultrasound Med. Biol., vol. 28, no. 1, pp. 101–113, 2002. [6] S. Langeland, J. Dapos, H. Torp, B. Bijnens, and P. Suetens, “A simulation study on the performance of different estimators for two dimensional velocity estimation,” in Proc. IEEE Ultrasonics Symp., vol. 2, 2002, p. 18591862. [7] R. Righetti, S. Srinivasan, and J. Ophir, “Lateral resolution in elastography,” Ultrasound Med. Biol., vol. 29, no. 5, pp. 695–7–4, 2003. [8] F. Viola and W. Walker, “A comparison of the performance of time delay estimators in medical ultrasound,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 50, p. 392401, 2003. [9] R. Zahiri and S. Salcudean, “Motion estimation in ultrasound images using time-domain cross-correlation with prior estimates,” IEEE Trans. Biomed. Eng., vol. 53, no. 10, pp. 1990–2000, 2006. 131 Chapter 5. High Frame Rate Ultrasound [10] G. Pinton, J. Dahl, and G. Trahey, “Rapid tracking of small displacements with ultrasound,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 53, p. 11031117, 2006. [11] R. Zahiri-Azar and S. E. Salcudean, “Time-delay estimation in ultrasound echo signals using individual sample tracking,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 55, no. 12, pp. 2640–265, 2008. [12] L. Taylor, B. Porter, D. Rubens, and K. Parker, “Three-dimensional sonoelastography: Principles and practices,” Phys. Med. Biol., vol. 45, pp. 1477–1494, 2000. [13] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulus imaging with 2-d transient elastography,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 49, no. 4, pp. 426–435, April 2002. [14] J. Bercoff, R. Sinkus, M. Tanter, and M. Fink, “Supersonic shear imaging: A new technique for soft tissue elasticity mapping,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 51, no. 4, pp. 396–409, 2004. [15] K. Nightingale, M. Palmeri, and G. Trahey, “Analysis of contrast in images generated with transient acoustic radiation force,” Ultrasound in Med. & Biol., vol. 32, no. 1, p. 6172, 2006. [16] J. J. Dahl, G. F. Pinton, M. L. Palmeri, V. Agrawal, K. R. Nightingale, and G. E. Trahey, “A parallel tracking method for acoustic radiation force impulse imaging,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 54, no. 2, pp. 301–312, 2007. [17] C. Schmitt, A. H. Henni, and G. Cloutier, “Characterization of timevarying mechanical viscoelastic parameters of mimicking deep vein thrombi with 2d dynamic elastography,” in Proc. of IEEE Ultrasonics Symposium, 2007, pp. 1009–1012. [18] A. H. Henni, C. Schmitt, and G. Cloutier, “Analytical modeling of plane shear wave diffraction by a radially layered cylinder for dynamic vascular elastography,” in Proc. of IEEE Ultrasonics Symposium, 2007, pp. 1713–1716. [19] S. Wang, W.-N. Lee, J. Luo, and E. E. Konofagou, “A composite imaging technique for high frame-rate and full-view cardiovascular ultrasound and elasticity imaging,” in Proc. of IEEE Ultrasonics Symposium, 2007, pp. 880–883. 132 Chapter 5. High Frame Rate Ultrasound [20] S. Wang, W.-N. Lee, J. Provost, J. Luo, and E. E. Konofagou, “A composite high-frame-rate system for clinical cardiovascular imaging,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 55, no. 10, pp. 2221–2233, 2008. [21] S. Brekke, C. B. Ingul., S. A. Aase., and H. G. Torp, “Increasing frame rate in ultrasound imaging by temporal morphing using tissue doppler,” in Proc. of IEEE Ultrasonics Symposium, 2004, pp. 118–121. [22] M. Fink, L. Sandrin, M. Tanter, S. Catheline, S. Chaffai, J. Bercoff, and J. Gennisson, “Ultra high speed imaging of elasticity,” in IEEE Ultrasonics Symposium, 2002, pp. 1811–1820. [23] C. Fabian, K. Ballu, J. Hossack, T. Blalock, and W. Walker, “Development of a parallel acquisition system for ultrasound research,” Proc. SPIE, vol. 4325, pp. 54–62, 2001. [24] T. Misaridis and J. A. Jensen, “Use of modulated excitation signals in medical ultrasound. part iii: High frame rate imaging,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 52, no. 2, pp. 208–219, 2005. [25] J. Shen and E. S. Ebbini, “Filter-based coded-excitation system for high-speed ultrasonic imaging,” IEEE Trans. Med. Imaging, vol. 17, no. 6, pp. 923–934, 1998. 133 Chapter 6 Conclusion 6.1 Summary The inverse problem of elastography was studied in this research [1–14]. The objective of the inverse problem is to find the elasticity of the tissue as a function of position, given the displacements as a function of position and time. In ultrasound elastography, usually only one component of the displacement, namely the axial component, is available from the measurements. Moreover, the measurement is restricted to the plane of imaging. In this case, the inverse problem is that of finding tissue elasticity based on partial displacement measurements. The wave equation approach to the inverse problem, in the case of partial displacement measurements, addresses the problem by finding the phase speed of the waves [2–5, 7, 8, 10, 11, 13, 14]. In this thesis, the limitations of this approach were clearly defined. In particular, it was shown that contrary to common assumptions, longitudinally polarized waves may not always travel with the so called shear wave speed in a medium. The shear wave speed is the phase speed of plane shear waves in an infinite medium and is given by µ cs = . (6.1) ρ It was shown that in a finite block of tissue, when a large profile mechanical exciter is used to compress the tissue longitudinally, the phase speed is given by µ cphase = 3 , (6.2) ρ which is the speed of longitudinal waves in thin rods. A further analysis of the phase speed revealed that the phase speed, in addition to the tissue elasticity, depends on the shape of the wave beam, and on the boundary conditions of the medium. Therefore a universal relationship between the phase speed and the tissue elasticity cannot be assumed. However in particular cases and under certain boundary conditions, the re134 Chapter 6. Conclusion lationship (6.1) holds and the wave speed can be used to estimate tissue elasticity. Different methods have been used for estimating the phase speed. These methods are based on local frequency estimation. For instance they measure the wavelength of the wave inside the medium. For accurate measurement of the wavelength, it is of interest to have a few periods of the wave observed in the length of the region of interest. To place more periods in the region of interest, the frequency of excitation should be increased. This is a common feature of all phase estimation methods; to achieve a higher resolution, a higher frequency of excitation should be used. In this thesis a technique was presented for increasing the frame rate of an ultrasound system to image higher frequencies of excitation. The conventional method of ultrasound imaging has a limited frame rate which depends on the depth of imaging, line density, and width of the region of interest. Typical frame rates are in the range of 50Hz to 200Hz. By the Nyquist criterion, this limits the frequency of excitation to 25Hz to 100Hz. The presented method in this thesis was easily implemented on a conventional ultrasound system, and increased the frame rate to a range of 1kHz to 10kHz. The presented high frame rate system was used to successfully track motions with a frequency of up to 500Hz. The tracked displacements were used to find the phase speed, and elasticity. The elastograms obtained from hard inclusion and soft inclusion phantoms clearly revealed the inclusion. At this point, it can be stated that all the primary objectives of this thesis were achieved. • A method was devised for accurately measuring the elasticity of tissue and tissue mimicking material samples for setting gold standards. • Tissue mimicking materials were fabricated with controllable elasticity. • The relationship between the phase speed and tissue elasticity was established through a theoretical analysis of the elastic wave equation. • An ultrasound elastography system was developed which uses the developed inversion algorithms for forming elastograms. Therefore the hypothesis that a wave equation based approached can be used in ultrasound elastography for forming elastograms has been confirmed. 135 Chapter 6. Conclusion 6.2 6.2.1 Contributions Rheometry A novel framework was devised in Chapter 2 for viscoelastic characterization of tissue and tissue mimicking material. Longitudinal wave excitation was used on blocks of test phantoms and displacements were imaged by an ultrasound system. Two methods were presented for analyzing the displacements. The first method estimated the Young’s modulus based on the distance from the peak to the node of the standing wave patterns inside the block. The second method used a model fitting scheme to find the values of the Young’s modulus and the relaxation time. The results were compared to conventional rheometry results obtained using force-displacement measurements. The error in the estimated elasticities is lowest for the range of 30-80kPa. The maximum error for the peak to node method was 11% and for the model fitting method was 20%. For the viscosities, only the trends could be compared as the measurements were over different frequency ranges. 6.2.2 Poly Vinyl Chloride Phantoms In Chapter 2 fabrication of homogeneous blocks of tissue mimicking material was reported from Super Soft PlasticTM , a Poly Vinyl Chloride based phantom using different proportions of the softener and hardener. Super Soft PlasticTM was shown to be a promising phantom material for elastography studies, as its Young’s (and hence shear) modulus can be adjusted over a wide range (15–150kPa), without changing its density and the speed of ultrasound cinf . Figure 6.1 summarizes the Young’s modulus obtained vs. the volume portion of hardener or softener used in the fabrication process for this material. However, the relaxation time and its frequency dependent behavior is not dependent on the volume portion of Softener or Hardener (refer to Fig. 2.5(b) in Chapter 2). Therefore without any additives, this material by itself is unsuitable for creating phantoms with significant viscosity contrast. 6.2.3 Phase Speed and Tissue Elasticity A theoretical analysis of the phase speed and its relationship to tissue elasticity was carried out in Chapter 3 which enhanced our understanding of the wave phenomena. From this analysis, it was concluded that the propagation of elastic waves trapped in finite media, such as soft tissue, is a complex phenomenon even without considering the viscous and nonlinear 136 Chapter 6. Conclusion Figure 6.1: Obtained Young’s modulus as a function of the volume portion of Plastic Softener or Hardener added during the fabrication of the phantoms effects. The term “wave” is generally applied to any vibration pattern inside the tissue. However the phase velocity cannot be defined for all the cases of the vibration patterns. Usually multiple waves traveling in different directions superimpose to create the vibration pattern. In this case the phase velocity might be well-defined for each wave, but not for the resultant of the superposition. Even in the simple case where a single frequency wave is traveling in a single direction, the geometry of the wave and the medium create multiple permissible phase speeds, and thus make it difficult to recover a meaningful phase speed for the traveling wave. Therefore no universal relationship can be proposed which would relate the phase speed to tissue elasticity. Based on the analysis carried out the following measures were recommended to minimize the artifacts in elastograms obtained using phase speed estimation. The main goal should be to create a single wave propagating in one direction at a single speed: • The excitation amplitude should be chosen small enough to reduce the reflected wave amplitudes to the level of other measurement noises. The damping present in the tissue helps by reducing the amplitude of the reflected waves from the boundaries of the tissue and bones. • The frequency of excitation should be kept to a minimum to prevent the appearance of higher modes. However lower frequencies result in lower damping, and therefore a trade off exists here. • The excitation pattern should be chosen so as to excite mainly the 137 Chapter 6. Conclusion lowest mode. 6.2.4 High Frame Rate Ultrasound System A novel high frame rate ultrasound system was introduced in Chapter 5 which could readily be implemented on conventional ultrasound systems without the need for any hardware modifications. The system uses the concept of sector subdivision for sequencing which is widely used in Doppler imaging to achieve high frame rates. For a sector of width one, i.e. a single scan-line, the frame rate is only limited by the time it takes for the deepest point to be scanned, i.e. limited by the speed of sound of 1540m/s in tissue. In addition, a novel phase correction scheme was presented to compensate for the time delays involved in the acquisition process. Necessary conditions on the frequency content of the tissue motion were presented in order for the displacements to be recovered properly. Since the devised system simply reorders the pulse sequence, the tissue is still imaged one line at a time. As a result, the total acquisition time with this system can be a fraction of a second. The probe and tissue should be stationary during this acquisition interval to prevent motion artifacts. The system can be used together with periodic excitations to study tissue characteristics, and as such is not meant as a replacement for other high speed imaging techniques which are used in applications with different requirements. 6.2.5 Ultrasound Elastography System The high frame rate ultrasound system was used to devise an ultrasound elastography system, which was reported in Chapter 5. The system uses periodic excitations, with multiple excitations frequencies to vibrate the tissue. The high frame rate ultrasound system images the tissue and provides the phase compensated phasor displacements at the frequencies of excitation. The phasor displacements are then used in an inversion algorithm, namely the algebraic inversion of the wave equation, to estimate the phase speed and hence tissue elasticity. The results obtained from different frequencies are averaged to give the final elastograms. The performance of the system was tested on hard inclusion and soft inclusion phantoms. The images obtained clearly show the geometry and location of the inclusions through changes in relative elasticity. However, measuring absolute values of elasticity has been a challenge which has not been fully addressed in this thesis. The underlying problem is that of estimating the local frequency (or wavelength) 138 Chapter 6. Conclusion of a function from a few samples of it at each of these samples. Using higher frequencies of excitation helps, but it also presents two new challenges. As the frequency goes higher, the damping of the wave becomes more severe in real tissue. Frequencies of above 500Hz cannot penetrate very deep inside real tissue. The second problem is the appearance of multiple modes and phase speeds. This problem could be addressed by using the curl of the displacements, as suggested in the next section. 6.3 Future Work In this thesis, a high-frame-rate elastography system was reported. The system works with periodic excitations in the range of 1-500Hz. The elasticity image update rate depends on the width of the ROI and the number of frames acquired and can be as low as a few images per second. The performance of the system was tested on tissue mimicking material. Along this line of research, a number of directions can be pursued: • Different excitation and inversion methods can be tested together with the high-frame-rate system to find an optimal combination for performing in vivo tests. • The parameters involved, such as the number and values of the frequencies of excitation, the frame rate of the high-frame-rate system, and the number of frames to be acquired could be optimized based on the application at hand. • The system could then be used in a clinical trial to prove its efficacy. Of particular interest are the breast and prostate elastography. The problem of local frequency estimation from a few samples of the signals is an inherently difficult problem. This is the problem of phase speed estimation for a small region of interest composed of for instance 100 pixels by 100 pixels. It is therefore of interest to circumvent the phase speed estimation, and use the wave equation directly to find the values of the mechanical properties. Along this line of research, a number of directions can be pursued: • The wave equation can be written for the shear component of the wave by application of the curl operator on the displacement data. To be able to use the curl operator, displacements should be measured in at least two directions and over a volume. The current research at the Robotics and Control Lab is addressing the problem of two directions. By using beam steering in conjunction with the high-frame-rate and phase compensation scheme of this thesis, lateral motions can be 139 Chapter 6. Conclusion measured accurately in addition to axial motions. • Acquisition of a 3D volume of beam steered ultrasound data with axial and lateral motion tracking using 3D ultrasound probes is a further step. • Different wave equation inversion algorithms based on the curl of the displacements can be implemented and their performance assessed on tissue mimicking material and real tissue. • Measuring absolute values of tissue elasticity is another challenging area which requires further investigation. 140 References [1] R. Sinkus, J. Lorenzen, D. S. amd M. Lorenzen, M. Dargatz, and D. Holz, “High-resolution tensor mr elastography for breast tumor detection,” Physics Med. Biol., vol. 45, pp. 1649–1664, 2000. [2] A. Manduca, T. Oliphant, M. Dresner, J. Mahowald, S. Kruse, E. Amromin, J. Felmlee, J. Greenleaf, and R. Ehman, “Magnetic resonance elastography: Non-invasive mapping of tissue elasticity,” Medical Image Analysis, vol. 5, pp. 237–2540, 2001. [3] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complexvalued stiffness reconstruction by for magnetic resonance elastography by algebraic inversion of the differential equation,” Magnetic Resonance in Medicine, vol. 45, pp. 299–310, 2001. [4] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulus imaging with 2-d transient elastography,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 49, no. 4, pp. 426–435, April 2002. [5] L. Sandrin, M. Tanter, J.-L. Gennisson, S. Catheline, and M. Fink, “Shear elasticity probe for soft tissues with 1-d transient elastography,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 49, no. 4, pp. 436– 446, April 2002. [6] J. Bercoff, S. Chaffai, M. Tanter, L. Sandrin, S. Catheline, M. Fink, J. L. Gennisson, and M. Meunier, “In vivo breast tumor detection using transient elastography,” Magnetic Resonance in Medicine, vol. 29, no. 10, pp. 1387–1396, 2003. [7] J. Bercoff, M. Muller, M. Tanter, and M. Fink, “Study of viscous and elastic properties of soft tissue using supersonic shear imaging,” in IEEE Ultrasonics Symposium, 2003, pp. 925–928. [8] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity and viscosity from measurement of shear wave speed dispersion,” J. Acoust. Soc. Am., vol. 115, no. 6, pp. 2781–2785, June 2004. 141 Chapter 6. Conclusion [9] J. Bercoff, R. Sinkus, M. Tanter, and M. FinkJ, “3d ultrasound-based dynamic and transient elastography: First in vitro results,” in IEEE Ultrasonics Symposium, 2004, pp. 28–31. [10] J. Bercoff, R. Sinkus, M. Tanter, and M. Fink, “Supersonic shear imaging: A new technique for soft tissue elasticity mapping,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 51, no. 4, pp. 396–409, 2004. [11] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelkaram, and J. Culiolic, “Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: An inverse problem approach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec 2004. [12] J. Fromageau, J.-L. Gennisson, C. Schmitt, R. L. Maurice, R. Mongrain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel mechanical properties with four ultrasound elastography methods and comparison with gold standard testing,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007. [13] J.-L. Gennisson, S. Catheline, S. Chaffa¨ı, and M. Fink, “Transient elastography in anisotropic medium: Application to the measurement of slow and fast shear wave speeds in muscles,” J. Acoust. Soc. Am., vol. 114, no. 1, pp. 536–541, July 2003. [14] A. Romano, J. Bucaro, P. Abraham, and S. Dey, “Inversion methods for the detection and localization of inclusions in structures utilizing dynamic surface displacements,” in Proceedings of SPIE Vol. 5503, vol. 5503, 2004, pp. 367–374. 142 Appendix A Conventional Rheometry The measurement of the viscoelastic properties of the tissue phantoms at the frequencies considered in this article (between 50Hz to 100Hz) by conventional rheometry (force-displacement measurements) is a nontrivial problem. Here we elaborate on this issue. In conventional compression based force-displacement rheometry, a block of test material with known dimensions w × l × h is placed between two parallel plates (see Fig. A.1). One of the plates is actuated and can move to compress the block. The force F needed to compress the material and the displacement of the moving plate d are measured. The stress and strain are found as, xx = τxx = d , h F . w×l (A.1) (A.2) The elasticity (stiffness) of the material is found as, E= τxx . (A.3) xx Conventional rheometry can also be used to study the behavior of the stiffness as a function of frequency. For this purpose the actuated plate of the device is excited by a time varying waveform, causing a time varying strain and a time varying stress. In this case, the stiffness would be S(jω) = F{τxx (t)} . F{ xx (t)} (A.4) where F() is the Fourier transform. The extension of static rheometry to a dynamic setting raises a number of very important issues which affect the measurement results. The first of these issues is the mounting of the force sensor. In the static case, the force sensor can be fitted between either of the plates and the rest of the device. In equilibrium the forces are different 143 Appendix A. Conventional Rheometry only by the weight of the test block and the weight of some parts of the force sensor. This difference can be measured and the device can be calibrated for this setting. In the dynamic case, however, there are significant issues related to each of the possible fittings. Assume that the force sensor is fitted between the actuator of the device and the moving plate (Fig. A.1). In this case the measured force by the sensor not only includes the applied force on the material block but also the inertial forces necessary to accelerate the plate and the mass of the block. The device needs to be calibrated, but the calibration is much more complicated than in the static setting. The inaccuracies involved in the measurements make this method very challenging to implement. Figure A.1: Simple schematic of the rheometry device in a conventional rheometry configuration. In this configuration the force sensor measures the force on the moving plate and moves with it. Another option is to fit the force sensor between the fixed plate and the chassis of the device (Fig. A.2). Although in this case the measured force is the force applied to the block to compress it from above, a second factor comes into play which complicates the problem. Due to the finite speed of propagation of longitudinal waves in the block, there would be a time delay (and hence a phase difference) between the measured force on the fixed end of the block and the force on the moving end of the block. The higher the frequency of the excitation, the larger the phase difference: φ(ω) = ωt = ω h . cfin (A.5) If the velocity of the longitudinal waves cfin were known, this formula could 144 Appendix A. Conventional Rheometry be used directly to compensate for the phase difference, S(jω) = F{τxx (t)} exp(jφ(ω)) , F{ xx (t)} (A.6) but this velocity is itself a function of the stiffness (2.30): cfin = E = ρ real(S(jω)) . ρ (A.7) This makes the measurement more challenging, but still feasible for low frequencies. Figure A.2: Simple schematic of the rheometry device in a conventional rheometry configuration. In this configuration the force sensor measures the force on the static plate and does not move. At higher frequencies, the problem becomes more severe. In fact equation (A.4) only applies to infinitesimal volumes of matter. Referring to Fig. A.3, this equation states that if the only forces acting on the infinitesimal volume are perpendicular forces to the top and bottom surfaces, the strain and stress in the vertical direction are proportional, the proportionality factor being the Young’s Modulus E. Whether this equation could be applied to a block of material depends on how similar the conditions are to the infinitesimal setting. As far as the frequency is concerned, in the infinitesimal setting, the height of the block is always negligible with respect to the wave length of a propagating wave. This is not always true for a block. In other words, the block cannot always be considered as a lumped element with lumped mechanical properties equal to that of its infinitesimal elements. If the wave 145 Appendix A. Conventional Rheometry Figure A.3: The infinitesimal volume which is used in the definition of Elastic modulus. length is comparable to the height of the block, equation (A.4) would no longer be valid for the block. One way to deal with the problem of the block not being a lumped element at higher frequencies is to reduce the height of the block, h. In other words, thin slabs could be used for rheometry measurements. But going back to the question of the similarity between the block and the infinitesimal settings, the boundary condition that “the only forces acting on the infinitesimal volume be perpendicular forces to the top and bottom surfaces” can never be achieved with the block; without any lubricant on the top and bottom plates of the device, the material would stick to the plates because of the friction and there would be tangential forces acting on the top and bottom surfaces of the block. The middle portion of the block, however, resembles the infinitesimal setting to a good extent. But reducing the block height h reduces the effective middle portion. Putting some lubricant helps but never removes the friction. It might be worth mentioning that this problem mainly appears for nearly incompressible material such as living tissue and tissue phantoms. These materials have a high tendency to expand laterally as they are compressed. From this discussion, it follows that rheometry of soft tissue by conventional compression type force-displacement measurement at higher frequencies is not a trivial problem. In fact the proposed ultrasound based rheometry technique can be considered as a simple solution for this problem. Most commercially available rheometry devices use shear rheometry but their bandwidth is limited to 50Hz or bellow. Even with higher bandwidth systems, similar issues would exist with soft tissue since the test blocks would not be lumped at higher frequencies. In this report we validated our estimated elasticity values with the values 146 Appendix A. Conventional Rheometry obtained from low frequency rheometry using the second configuration (Fig. A.2). The relaxation times however could not be validated directly, except for the general trends. 147 Appendix B Tissue Mimicking Materials for Ultrasound Elastography In this appendix, a brief review of some of the most commonly used tissue mimicking materials is presented. Also a detailed recipe for fabricating PVC phantoms is given. B.1 Literature Review A comprehensive study of hydrogels, namely gelatine and agar based phantoms has been carried out in [1]. The purpose of the study is assessment of the acoustic and mechanical properties of these phantoms for ultrasound elastography. The authors mainly study Young’s modulus and sound speed. They give exact recipes for preparing the phantoms, and carry out different measurements and report the results. The phantoms are based on gelatine or agar solution in a mixture of dionized water and n-propanol. n-propanol is used to increase the sound speed. Graphite or glass micro-beads are used as scatterers. Some of the methods they use in preparation of the phantoms, such as use of vacuum chambers for degasing the mixtures, or rotating the mixtures to keep the graphite afloat while the gel congeals are not really appropriate for making complex anatomical phantoms, as they might be for cylinderical or cubic phantoms. The authors also study the effect of time passage on the degradation of the prepared phantoms. Another study on phantoms appropriate for ultrasonic elastography is carried out by Madsen et.al. [2]. The authors perform an informative literature survey on the state of the art phantoms used in ultrasonography and elastography. They note that most of the these studies lack a tracking of the aging of the phantoms and the changes in their characteristics. As a result the phantoms are not stable during the development of imaging techniques and hence these researches suffer a methodical problem. They also propose gelatine based phantoms with dispersed safflower oil for ultrasonic elastography. Exact recipes are given. Different proportions of ingredients are used 148 Appendix B. Tissue Mimicking Materials to achieve different characteristics. The group studies the stability of acoustic and mechanical properties of the phantoms, over an extended period of time in [3]. They also study the stability issue for agar and gelatine based phantoms without oil dispersions [4]. A comprehensive literature survey has been carried out by Hassan and Peppas on Poly(vinyl alcohol) polymers in [5]. This material has a wide variety of applications in biomedical engineering from contact lenses to artificial implants, such as artificial cartilages. Therefore numerous methods have been developed for specific applications. The paper does not aim to give a certain recipe for preparation of Poly(vinyl alcohol) gels, rather gives an overview of different methods and material used for obtaining different gels. Another seminal paper has recently been published by Madsen et.al. [6]. In this paper the authors present the fabrication of two anthropomorphic breast phantoms for Ultrasound and MRI Elastography. Six tissue types have been prepared. The material used is gelatine with dispersed safflower oil. A well-detailed recipe is given as well as the involved molding procedure used to prepare the phantom. The paper sets a good standard for preparing anthtopomorphic phantoms. The only downside of the fabricated phantom is that it should be kept in safflower oil, even during experiments to prevent desiccation. In [7] the authors present a bone-mimicking material base on plaster, CMC vegetable adhesive and chopped flax which shows attenuation and speed characteristics similar to bovine bone. Another bone mimicking tissue has been proposed in [8] using apatite and PVA. Sound speed, characteristic impedance and the density of the phantom matches those of the bovine cancellous bone. In [9] the authors use a 13% gelatine, 10% graphite, 10% formaldehyde phantom with an embedded cellophane locator marker and vibrate it to create strain-rate (SR) images using doppler-ultrasound. The phantom is based on the work of Hall et al. in [1]. The claim is that the mixture closely simulates the stiffness and scattering properties of soft tissue. The authors of [10] investigate the temperature and frequency dependence of four types of commercially available phantoms: urethane rubber from ATS Labs (Bridgeport, CT), two types of condensed milk from Gammex-RMI (Middleton, WI) and ZerdineTM from CIRS Inc. (Norfolk, VA). They also study an agar based phantom developed through a European Commission (EC)-funded project. They mention that the International Electrotechnical Commission 1390 (IEC 1996) and American Institute of Ultrasound in Medicine Standard 1990 (AIUM 1990) standards for TMMs 149 Appendix B. Tissue Mimicking Materials recommend an acoustic velocity of 1540 m/s, attenuation coefficients of 0.5 dB/(cm MHz) and 0.7 dB/(cm MHz) for the frequency range 2 to 15 MHz with a linear response of attenuation to frequency, f 1 . PolyVinyl Alcohol Cryogel (PVA-C) phantoms which have undergone different numbers of freeze and thaw cycles to adjust their acoustic and mechanical properties are used for vessel elastography in [11]. The authors cite [12] as their reference for building the phantoms, and do not mentioned the exact method they have used, except the number of freeze and thaw cycles. In [13] PVA based phantoms are proposed which have optical, acoustic and mechanical properties similar to breast tissue. The authors are performing Ultrasound Assisted Optical Elastography (UAOE), a new approach to elastography. They describe in detail the recipe for making the PVA phantoms and the freeze-n-thaw cycles used. Their work is based on [5] and [14]. They also describe methods for measuring optical, acoustic and mechanical properties of the phantoms. The measured properties for two different PVAs and 5 different numbers of freeze-n-thaw cycles are given. Simple PVA based material are used to fabricate a number of antropomorphic phantoms in [15]. They have used usuall freeze-n-thaw cycles to increase the stiffness of their gels. They also study the effect of adding color to dye the phantom on some of the characteristics. Complex tissue mimicking phantoms are proposed in [16] for multimodal imaging of prostate. These contain water, agarose, lipid particles, protein, Cu++, EDTA, glass beads, and thimerosal (preservative). Nevertheless the mechanical properties of the phantoms are not of concern to the authors. Nonlinear mechanical properties of an agar-based and a gelatine-based phantom under large static strains are reported in [17]. The authors used bleach for preservation and graphite to create scattering and attenuation. Exact percentages used are mentioned but not the recipe. In [18] five different gelatine-based phantoms are developed for assesment of tissue elasticity and viscosity within the framework of shear wave elasticity imaging and acoustic remote palpation. Glycerol was added to gelatine to increase the viscosity of the phantoms. Shear wave speed, speed of sound and attenuation coefficient for a 1MHz probe are given for the phantoms. Quantity of ingredients are given but not the recipe. The fabrication of an anthropomorphic phantom of the oesophagus is reported in [19]. The work is of particular interest, as it describes the methods used for integration of the anatomical features as well as the recipe for preparing them. Different tissue-mimicking parts are fabricated using different consentrasions of Al2 O3 and SiC particles inside an agar-based hydrogel. 150 Appendix B. Tissue Mimicking Materials A mixture of water, glycerol and benzalkonium is used to dissolve the agar. Exact quantity of the material used and the recipe are given. Another hydrogel which has been proposed for making ultrasonic elastography phantoms is polyacrylamide [20]. The authors have used TiO2 to adjust the acoustic properties of the phantom, and different concentrations of acrylamide to alter the elasticity. The recipe for using this material is given. The negative point about this material is that although polyacrylamide is nontoxic, but since the polymerization is not perfect, some residual acrylamide is always present, and this latter is a neurotoxin. Nevertheless it seems that small quantities of acrylamide is present in fried food and wether it is a carcinogen is still a point of debate. Special handling practices are recommended. Another PVA based phantom for photoacoustic mammography has been developed in [14]. In this imaging technique, the ultrasonic waves are produced internally by local heating of the tissue with near-infrared laser pulses. In a first method, the authors use freeze-n-thaw cycles to achive at the desired photoacoustic properties. In a second method they use a mixture of water and dimethylsulphoxide to dissolve the PVA, and reach at a transparent gel. Then by addition of borosilicate glass microspheres they alter the turbidity of the phantom. Except for the densities no study of mechanical properties is carried out. Kondo and Fujimoto have reported oil-gel phantoms [21]. The main advantage of these phantoms is that they do not rust. The phantoms are based on ethylene glycol (car radiator antifrost) and propylene glycol or polypropylene glycol. Etylene glycol is toxic in that it must not be ingested. Propylene glycol is safe and approved by FDA. It is used in foods and sexual lubricants. The author only report sound speed, density and some relative measure of hardness. No recipes are given. Madsen et.al. propose a phantom geometry for testing Ultrasound and NMR elastography methods in [22]. The phantom consists of a series of spherical balls of different sizes embedded in a background. Three different phantoms are fabricated with this geometry. Two of them are agar-gelatine based and the third one involves dispersions of different concentrations of microscopic oil droplets in a gelatin matrix. The exact percentages of the material are given, but for the recipe, the authors refer the reader to [2] and [4]. The paper mainly concentrates on the molding procedure used to obtain the phantoms. Polyvinyl chloride(PVC) is another polymer which has the potential of being used as a base for elastography and ultrasonography phantoms. To the best of our knowledge, this material has not been investigated thoroughly 151 Appendix B. Tissue Mimicking Materials in the literature. Although some reports on the use of this material exist [23] [24] [25] and [26]. B.2 B.2.1 PVC Phantom Recipe Introduction Super Soft PlasticTM phantoms are made by mixing appropriate proportions of Super Soft PlasticTM with either Plastic SoftenerTM or PlasticTM . The ultrasound speed is pretty much independent of the mixture ratio and is around 1420m/s, slightly lower than the ultrasound speed in real tissue which is 1540m/s. The amount of softener or hardener however determines the elasticity or stiffness of the final phantom. Without adding any particles, the phantoms would be transparent to ultrasound. Cellulose with average 50micron particle size is added to create scatterer pattern. The ratio of scatterer to plastic volume should not exceed 4%, otherwise excessive scattering would limit the penetration depth of the ultrasound to a few centimeters. Good results are obtained using 1% cellulose which is about 1/4tsp for 100mL plastic. B.2.2 Preparation Steps The following steps should be taken. • The plastic, softener and hardener container should be shaken well before pouring the plastic, since the plasticizer and the liquid plastic tend to form a two-phase mixture. • The plastic, softener or hardener is measured and poured inside a beaker. • The beaker is heated over the hotplate while trying to maintain a uniform temperature distribution in the liquid by mixing. • When the plastic turns into a transparent liquid, the cellulose is added slowly while stirring to prevent creation of lumps. • Let the liquid stay at a high enough temperature to keep it in the liquid state to allow for the air bubbles to surface. This may take 10 to 30 minutes. • Pour the liquid slowly into the mold to allow for the generated air bubbles during the pouring process to remain on the surface and not get buried. • Let the mold to cool down well before trying to remove the phantom. Remove the phantom slowly by trying to disconnect its surfaces from 152 Appendix B. Tissue Mimicking Materials the mold surfaces. B.2.3 Other Important Notes and Tips The following notes have been made from several attempts at making PVC phantoms: • Since the liquid plastic is heated close to its evaporation point, a lot of fume is generated which may pose health and safety risks to you and the other people working in the room. The whole process of heating, molding and cooling of the phantoms should therefore be performed under a fume hood. The company providing the material sells them for making fish lures. They recommend fabrication under the kitchen fume hood. I would recommend using a chemical lab fume hood for complete protection. • There is no water-based solvent for PVC. In other words commercial soaps, washing liquids and powders cannot clean the lab equipment, your hands, or anything which comes in touch with the PVC plastics. Acetone (finger polish remover) is a chemical that can be used. Acetone is very evanescent however and its vapor also poses health and safety risks. The surface to be cleaned should be wiped with acetone. • Hotplates come in a number of different types. The ones with controllable plate temperature and an even temperature distribution tend to be the best. However since these hotplates are not available in every case, what happens is that the plate temperature usually goes above the necessary temperature for heating the plastic. If the beaker is directly placed over the hotplate, this may cause the liquid at the bottom of beaker to burn and stick to the beaker bottom and also darkens the color of the phantom. To get a milder temperature distribution in this case a pan can be placed between the plate and the beaker. A better method would be filling a pot with sand and dipping the beaker inside the sand and heating the pot. This would create a more even temperature distribution inside the beaker. • When retaining the plastic in the molten state to let the air bubbles to surface, you don’t want the temperature to go very high, because the lighter chemicals would evaporate and change the properties of your phantom. Also the phantom would start to darken too much and loose its visual appeal. On the other hand too low a temperature would create a very viscose environment in which small air bubbles need ages to surface. So a compromise should be achieved by experience and looking at how fast the air bubbles are surfacing. 153 Appendix B. Tissue Mimicking Materials • The temperature at which the molten plastic is kept to allow for the air bubbles to surface depends on the ratio of hardener or softener to plastic. Generally a harder plastic requires a higher temperature to be liquid enough for the air bubbles to surface than a softer plastic. The temperature range is from 150C to 180C. So 150C is used for a very soft plastic and 180C for a very hard one. • When the liquid touches the cold metal surfaces of a mold, it will immediately cool down and coagulate. The best practice would be to heat up the mold to the same temperature as liquid before pouring in the molten plastic. Otherwise, care should be taken in pouring the molten plastic. Too slow a rate would create a laminated phantom because of the cooling down of the previously poured material. On the other had too fast a rate would trap the air bubbles in the phantom. • The ratio of softener or hardener to plastic volume should not exceed one because the created phantom would be too soft or too hard and very difficult to handle. 2/3 limit would be a more conservative option creating a large elasticity interval (about one decade). • During the heating process of the liquid plastic, at the gel point solid masses start to form inside the solution which would impede the stirring process. If the magnetic stirrer is used a lower stirring rate could be used until the solid masses are melted and a uniform liquid is achieved. Generally speaking the stirring should be kept to a minimum as it sucks in air bubbles. However without stirring it would be difficult to maintain an even temperature distribution in the liquid. 154 References [1] 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, 1997. [2] E. Madsen, G. Frank, T. 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. [3] E. L. Madsen, M. A. Hobson, H. Shi, T. Varghese, and G. R. Frank, “Stability of heterogeneous elastography phantoms made from oil dispersions in aqueous gels,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 261–270, 2006. [4] ——, “Tissue-mimicking agar/gelatin materials for use in heterogeneous elastography phantoms,” Physics in Medicine and Biology, vol. 50, p. 55975618, 2005. [5] C. M. Hassan and N. A. Peppas, “Structure and applications of poly(vinyl alcohol) hydrogels produced by conventional crosslinking or by freezing/thawing methods,” Advances in Polymer Science, vol. 153, pp. 37–65, 2000. [6] E. L. Madsen, M. A. Hobson, G. R. Frank, H. Shi, J. Jiang, T. J. Hall, T. Varghese, M. M. Doyley, and J. B. Weaver, “Anthropomorphic breast phantoms for testing elastography systems,” Ultrasound in Medicine and Biology, vol. 32, no. 6, pp. 857–874, 2006. [7] H. Asaba, E. Ohdaira, N. Masuzawa, and M. Ide, “Fundumental study to develop bone-mimicking phantom,” Japanese Journal of Applied Physics, vol. 38, no. 5B, pp. 3412–3413, May 1999. [8] T. Hirai, E. Ohdira, N. Masuzawa, K. Itoh, and T. Matozaki, “Fundamental study of bone-mimicking phantom using apetite,” Japanese Journal of Applied Physics, vol. 40, no. 5B, pp. 3905–3906, May 2001. 155 Appendix B. Tissue Mimicking Materials [9] M. Belohlavek, V. B. Bartleson, and M. E. Zobitz, “Real-time strain rate imaging:validation of peak compression and expansion rates by a tissue-mimicking phantom,” Echocardiography: A Jrnl. of CV Ultrasound and Allied Tech, vol. 18, no. 7, pp. 565–571, October 2001. [10] J. E. Browne, K. V. Ramnarine, A. J. Watson, and P. R. Hoskins, “Assessment of the acoustic properties of common tissue-mimicking test phantoms,” Ultrasound in Med. & Biol., vol. 29, no. 7, pp. 1053–1060, 2003. [11] E. Brusseau, P. Delachartre, and D. Vray, “Axial strain imaging of vessel mimicking cryogel phantoms,” in Proceedings of Ultrasonics Symposium 2000 IEEE, vol. 2. IEEE, October 2000, pp. 1817–1820. [12] K. C. Chu and B. K. Rutt, “Polyvinyl alcohol cryogel: An ideal phantom material for mr studies of arterial flow and elasticity,” Magnetic Resonance in Medicine, vol. 37, no. 2, pp. 314–319, Feb 1997. [13] C. U. Devi, R. M. Vasu, and A. K. Sood, “Design, fabrication, and characterization of a tissue-equivalent phantom for optical elastography,” Journal of Biomedical Optics, vol. 10, no. 4, p. 044020, July-August 2005. [14] A. Kharine, S. Manohar, R. Seeton, R. G. Kolkman, R. A. Bolt, W. Steenbergen, and F. F. de Mul, “Poly(vinyl alcohol) gels for use as tissue phantoms in photoacoustic mammography,” Physics in Medicine and Biology, vol. 48, pp. 357–370, 2003. [15] K. Surry, H. Austin, A. Fenster, and Peters, “Poly(vinyl alcohol) cryogel phantoms for use in ultrasound and mr imaging,” Physics in Medicine and Biology, vol. 49, p. 55295546, 2004. [16] W. D. D’Souza, E. L. Madsen, O. Unal, K. K. Vigen, G. R. Frank, and B. R. Thomadsen, “Tissue mimicking materials for a multi-imaging modality prostate phantom,” Med. Phys., vol. 28, no. 4, pp. 688–700, April 2001. [17] R. Q. Erkamp, A. R. Skovoroda, S. Y. Emelianov, and M. O’Donnell, “Measuring the nonlinear elastic properties of tissue-like phantoms,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 51, no. 24, pp. 410–419, April 2004. 156 Appendix B. Tissue Mimicking Materials [18] S. Girnyk, A. Barannik, E. Barannik, V. Tovostiak, 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 Med. and Biol., vol. 32, no. 2, p. 211219, 2006. [19] S. Inglis, K. V. Ramnarine, J. N. Plevris, and W. N. McDicken, “An anthropomorphic tissue-mimicking phantom of the oesophagus for endoscopic ultrasound,” Ultrasound in Med. and Bio., vol. 32, no. 2, pp. 249–259, 2006. [20] K. ichi Kawabata, Y. Waki, T. Matsumura, and S. ichiro Umemura, “Tissue mimicking phantom for ultrasonic elastography with finely adjustable elastic and echographic properties,” in Proceedings of 2004 IEEE International Ultrasonics, Ferroelectrics, and Frequency Control Joint 50th Anniversary Conference. IEEE, 2004, pp. 1502–1505. [21] T. Kondo and H. Fujimoto, “Ultrasound tissue-mimicking materials using oil gel and measurement of their characteristics,” Japanese Jounral of Applied Physics, vol. 41, no. 5B, pp. 3598–3599, 2002. [22] E. L. Madsen, G. R. Frank, M. A. Hobson, H. Shi, J. Jiang, T. Varghese, and T. J. Hall, “Spherical lesion phantoms for testing the performance of elastography systems,” Physics in Medicine and Biology, vol. 5, pp. 5983–5995, 2005. [23] B. Chan, J. V. Merton-Gaythrope, M. P. Kadaba, S. Zafaranloo, and D. Bryk, “Acoustic properties of polyvinyl chloride gelatin for use in ultrasonography,” Radiology, pp. 215–216, 1984. [24] E. J. Chen, J. Novakofski, K. Jenkins, and W. D. O. Jr., “Ultrasound elasticity measurements of beef muscle,” in Proc. of IEEE 1994 Ultrasonics Symposium, 1994, pp. 1459–1462. [25] E. J. Chen, J. Novakofski, W. K. Jenkins, and W. D. O. Jr., “Youngs modulus measurements of soft tissues with application to elasticity imaging,” IEEE Transcations on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 43, no. 1, pp. 191–194, Jan. 1996. [26] T. Podder, D. Clark, J. Sherman, D. Fuller, E. Messing, D. Rubens, J. Strang, Y. Zhang, W. O’Dell, W. Ng, and Y. Yu, “Effects of tip geometry of surgical needles: An assessment of force and deflection,” in Proc. of the 3rd European Medical and Biological Engineering Conference, Prague, Czech Republic, Nov. 2005. 157 Appendix C Phasor Calculations Consider a scenario in which the tissue is continuously excited by a vibrating mechanical exciter. Practical limitations such as the sampling frequency, the finite time of observation, and the synchronization of the imaging and excitation systems impose constraints on the (time) shape of the excitation. The excitation is usually chosen to be a periodic signal. Without a loss of generality, assume that the excitation is sinusoidal at a frequency of fe . When multiple frequencies exist in the excitation, the exact same analysis applies for each of the components if linearity is assumed. In the steady state, the displacement of the point x ¯ = (x, y, z) in the tissue at time t can be written as, ui (¯ x, t) = Ai (¯ x) exp(j2πfe t) exp(jφi (¯ x)) +Ai (¯ x) exp(−j2πfe t) exp(−jφi (¯ x)) (C.1) i = x, y, z where u ¯ = (ux , uy , uz ) is the displacement vector. Since the two terms on the right hand side are complex conjugates, it is sufficient to use one for the ¯ = (Ux , Uy , Uz ) defined by analysis. We call U Ui (¯ x) = Ai (¯ x) exp(jφi (¯ x)) i = x, y, z (C.2) the phasor profile of the displacements at the frequency of excitation fe . Note that the phasor profile does not have a time dependency. To compute it, however, the displacements must be observed over an interval of time. ¯ can be Denoting the Fourier transform in the time variable by Ft {·}, U computed from u ¯, ¯ (¯ U x) = Ft {¯ u(¯ x, t)}|ω=2πfe (C.3) The phasor profile is used in the inversion techniques to estimate the tissue elasticity. A group of the techniques are based on the direct inversion of the wave equation. In some of these techniques, it is assumed that a component of the displacements satisfies the shear wave equation. For instance 158 Appendix C. Phasor Calculations for the y-component, one has ∂ uy (¯ x, t) = c2s (¯ x)∇2 uy (¯ x, t) ∂t2 (C.4) x)/ρ is the classically known shear wave speed and ∇2 where cs (¯ x) = µ(¯ is the Laplacian, ∂2 ∂2 ∂2 ∇2 = + + (C.5) ∂x2 ∂y 2 ∂z 2 The equation for the phasor profile can be found by taking the Fourier transform of (C.4) in the time variable, and evaluating it at the frequency of excitation as in (C.3), −(2πfe )2 Uy (¯ x) = cs (¯ x)2 ∇2 Uy (¯ x) (C.6) The shear modulus can be computed from µ(¯ x) = ρc2s (¯ x) = −ρ(2πfe )2 Uy (¯ x) 2 ∇ Uy (¯ x) (C.7) Another technique which is used for the inversion is the phase gradient technique. Assume that a plane wave is propagating in the n ˆ direction in the medium. Assume that the coordinate axes have been chosen, so that the direction of particle motion for the wave is in the y-direction. The y-component of the displacement can be written as 1 n ˆ·x ¯)) λ = A exp(j2π(fe t)) exp(−jkˆ n·x ¯) uy (¯ x, t) = A exp(j2π(fe t − (C.8) where λ is the wave length and k is the wave number, k= 2π 2πfe = λ cs (C.9) Note that for brevity the conjugate term present in (C.2) has been dropped in (C.8). The phasor profile for this traveling wave becomes, Uy (¯ x) = A exp(−jkˆ n·x ¯) (C.10) 159 Appendix C. Phasor Calculations The phase of the phasor profile is given by, ∠Uy (¯ x) = −kˆ n·x ¯ + ∠A (C.11) To calculate the shear modulus from the phasor profile, the wave number k is first found. This is done by taking the gradient of the phase in the n ˆ direction, k = −ˆ n · ∇∠Uy (¯ x) (C.12) The shear wave speed is then given by, cs = − 2πfe n ˆ · ∇∠Uy (¯ x) (C.13) and the shear modulus by, µ=ρ 2πfe n ˆ · ∇∠Uy (¯ x) 2 (C.14) 160
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A wave equation approach to ultrasound elastography
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A wave equation approach to ultrasound elastography Baghani, Ali 2010
pdf
Page Metadata
Item Metadata
Title | A wave equation approach to ultrasound elastography |
Creator |
Baghani, Ali |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | Tissue elastography is a field of medical imaging which deals with obtaining images of tissue mechanical properties, in particular tissue elasticity. Every elastography imaging system has three major components: an exciter, an imaging system, and an image processing software module. The principle of operation of these systems is that softer and harder tissues react differently to a mechanical excitation. The imaging system captures images of the tissue as the mechanical exciter deforms the tissue. The image processing software module then computes the tissue motion from the acquired images, and solves an inverse problem to obtain tissue elasticity from tissue displacements. In ultrasound elastography all three components of the displacement are not available, and the measurement is carried out only on the imaging plane. The problem of inversion using partial measurements is an inherently challenging problem. To solve this problem the phase speed is estimated in the first place. This thesis studies the process of estimating tissue elasticity from tissue displacements, using the phase speed estimation. It is shown that, in general, the elasticity cannot be inferred from the phase speed, but under certain boundary and excitation conditions, a relationship exists between the phase speed and the elasticity. However the relationship depends on the boundary and excitation conditions. Based on the developed results, a novel rheometry method is proposed for visco-elastic characterization of tissue samples. To increase the resolution of the elasticity images it is of interest to use higher frequencies of excitation. However the inherent low frame rate of ultrasound systems posed a limitation. A novel high-frame-rate ultrasound system is introduced which is capable of tracking motions of up to 500Hz. The high-frame-rate system is used in conjunction with the inversion algorithms to form an elastography system. The performance of the system is tested experimentally on tissue mimicking material having hard and soft inclusions. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-12 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0070904 |
URI | http://hdl.handle.net/2429/18032 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2010-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2010_spring_baghani_ali.pdf [ 10.86MB ]
- Metadata
- JSON: 24-1.0070904.json
- JSON-LD: 24-1.0070904-ld.json
- RDF/XML (Pretty): 24-1.0070904-rdf.xml
- RDF/JSON: 24-1.0070904-rdf.json
- Turtle: 24-1.0070904-turtle.txt
- N-Triples: 24-1.0070904-rdf-ntriples.txt
- Original Record: 24-1.0070904-source.json
- Full Text
- 24-1.0070904-fulltext.txt
- Citation
- 24-1.0070904.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0070904/manifest