A Wave Equation Approach toUltrasound ElastographybyAli BaghaniB.A.Sc., Sharif University of Technology, 2004M.A.Sc., University of Tehran, 2006A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinThe Faculty of Graduate Studies(Electrical and Computer Engineering)The University of British Columbia(Vancouver)January, 2010c© Ali Baghani 2010AbstractTissue elastography is a field of medical imaging which deals with obtainingimages of tissue mechanical properties, in particular tissue elasticity. Everyelastography imaging system has three major components: an exciter, animaging system, and an image processing software module. The principleof operation of these systems is that softer and harder tissues react differ-ently to a mechanical excitation. The imaging system captures images of thetissue as the mechanical exciter deforms the tissue. The image processingsoftware module then computes the tissue motion from the acquired images,and solves an inverse problem to obtain tissue elasticity from tissue displace-ments. In ultrasound elastography all three components of the displacementare not available, and the measurement is carried out only on the imagingplane. The problem of inversion using partial measurements is an inherentlychallenging problem. To solve this problem the phase speed is estimated inthe first place. This thesis studies the process of estimating tissue elasticityfrom tissue displacements, using the phase speed estimation. It is shownthat, in general, the elasticity cannot be inferred from the phase speed, butunder certain boundary and excitation conditions, a relationship exists be-tween the phase speed and the elasticity. However the relationship dependson the boundary and excitation conditions. Based on the developed results,a novel rheometry method is proposed for visco-elastic characterization oftissue samples. To increase the resolution of the elasticity images it is ofinterest to use higher frequencies of excitation. However the inherent lowframe rate of ultrasound systems posed a limitation. A novel high-frame-rate ultrasound system is introduced which is capable of tracking motionsof up to 500Hz. The high-frame-rate system is used in conjunction withthe inversion algorithms to form an elastography system. The performanceof the system is tested experimentally on tissue mimicking material havinghard and soft inclusions.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiStatement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.3 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . 7References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Measurement of Viscoelastic Properties of Tissue Mimick-ing Material Using Longitudinal Wave Excitation . . . . . 152.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2 Theory of Wave Motion . . . . . . . . . . . . . . . . . . . . . 172.2.1 Longitudinal Plane Waves in an Infinite Medium . . . 192.2.2 Longitudinal Plane Waves in a Finite Medium . . . . 202.3 Experimental Methods and Setup . . . . . . . . . . . . . . . 222.4 Proposed Rheometry Methods . . . . . . . . . . . . . . . . . 252.4.1 Peak to Node Method . . . . . . . . . . . . . . . . . . 282.4.2 Model Fitting Method . . . . . . . . . . . . . . . . . 292.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30iiiTable of Contents2.5.1 Peak to Node Results . . . . . . . . . . . . . . . . . . 302.5.2 Model Fitting Results . . . . . . . . . . . . . . . . . . 302.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.7 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . 39References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413 Theoretical Limitations of the Elastic Wave Equation Inver-sion for Tissue Elastography . . . . . . . . . . . . . . . . . . . 463.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.2 Local Inversion of the Wave Equation . . . . . . . . . . . . . 493.3 A Critique on the Terminology Associated with the WaveEquation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.4 Phase Speed of Mechanical Waves in Elastic Mediums . . . . 553.4.1 Dependence of the Phase Speed on the Geometry ofthe Wave: Infinite Mediums . . . . . . . . . . . . . . 553.4.2 Dependence of the Phase Speed on the Geometry ofthe Medium: Wave Guides . . . . . . . . . . . . . . . 593.5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . 613.5.1 Shear Waves in a Tissue Mimicking Infinite Medium . 613.5.2 Waves in a Tissue Mimicking Cylindrical Rod . . . . 643.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704 The Influence of the Boundary Conditionson Longitudinal Wave Propagationin a Viscoelastic Medium . . . . . . . . . . . . . . . . . . . . . 784.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.2 Governing Equations in a Linear Viscoelastic Medium . . . . 804.3 Wave Equation and the Propagation Speed . . . . . . . . . . 824.3.1 Longitudinal Wave in a Finite Medium Approximat-ing a Thin Rod . . . . . . . . . . . . . . . . . . . . . 834.3.2 Longitudinal Wave Induced by a Point Source Exciter 844.4 Determining the Wave Speed from the Spectra of the Dis-placements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.4.1 Case (i) . . . . . . . . . . . . . . . . . . . . . . . . . . 894.4.2 Case (ii) . . . . . . . . . . . . . . . . . . . . . . . . . 904.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 914.5.1 Dynamic Finite Element Analysis . . . . . . . . . . . 91ivTable of Contents4.5.2 Simulation of the Wave Propagation . . . . . . . . . . 934.5.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . 964.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 984.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045 A High Frame Rate Ultrasound System for the Study ofTissue Motions . . . . . . . . . . . . . . . . . . . . . . . . . . . 1075.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 1075.2 Sector-Based High Frame Rate Acquisition Technique . . . . 1095.2.1 Sequencing . . . . . . . . . . . . . . . . . . . . . . . . 1095.2.2 Time Delay Compensation . . . . . . . . . . . . . . . 1135.2.3 Implementation Issues . . . . . . . . . . . . . . . . . . 1165.3 Requirements on the Tissue Motion . . . . . . . . . . . . . . 1195.4 Applications and Experimental Results . . . . . . . . . . . . 1225.5 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . 128References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1316 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1346.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1346.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 1366.2.1 Rheometry . . . . . . . . . . . . . . . . . . . . . . . . 1366.2.2 Poly Vinyl Chloride Phantoms . . . . . . . . . . . . . 1366.2.3 Phase Speed and Tissue Elasticity . . . . . . . . . . . 1366.2.4 High Frame Rate Ultrasound System . . . . . . . . . 1386.2.5 Ultrasound Elastography System . . . . . . . . . . . . 1386.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 139References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141AppendicesA Conventional Rheometry . . . . . . . . . . . . . . . . . . . . . 143B Tissue Mimicking Materials for Ultrasound Elastography 148B.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . 148B.2 PVC Phantom Recipe . . . . . . . . . . . . . . . . . . . . . . 152B.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 152vTable of ContentsB.2.2 Preparation Steps . . . . . . . . . . . . . . . . . . . . 152B.2.3 Other Important Notes and Tips . . . . . . . . . . . . 153References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155C Phasor Calculations . . . . . . . . . . . . . . . . . . . . . . . . 158viList of Tables2.1 Fabrication and measured properties of the five test phantoms 242.2 Parameters of the ultrasound imaging system . . . . . . . . . 242.3 Parameters of the speckle tracking motion estimation algorithm 252.4 Estimated values of longitudinal wave speed cfin and Young’smodulus for the test phantoms using the peak to node method 322.5 Estimated values of longitudinal wave speed cfin and Young’smodulus for the test phantoms using the model fitting method 332.6 Summary of the elastic properties of the five test phantoms . 373.1 Mechanical property values used for the simulation . . . . . . 613.2 Mechanical property values used for the simulation . . . . . . 644.1 The theoretical and estimated wave speed for three differentconditions. The estimated wave speed are based on threemethod: first, the frequency that the peak of the spectrumhappens was used to reconstruct the wave speed; second, theprevious value was enhanced to account for the effect of vis-cosity; third, the wave-front due to a transient step was ana-lyzed to recover the speed. . . . . . . . . . . . . . . . . . . . . 974.2 The theoretical and estimated wave speed for an experimentwith a soft PVC phantom. A flat and a point-source excita-tion were applied to explore the effect of boundary conditionson the speed of longitudinal waves. . . . . . . . . . . . . . . . 985.1 Parameters of the high frame rate system . . . . . . . . . . . 124viiList of Figures2.1 Simple schematic of the rheometry device . . . . . . . . . . . 232.2 A typical normalized displacement phasor profile, belongingto phantom No.3 at a frequency of 90Hz. The (a) magnitudeand (b) phase of the phasor profile (c) the same normalizedphasor profile in the complex plane. The profile is normalizedto the phasor of the point at x = 40mm. . . . . . . . . . . . . 262.3 Fitted model to the measured displacement data. The data isthe same as that of Fig. 2.2, which is from phantom No. 3 at90Hz. (a) absolute value and (b) phase of the phasor profile(c) the same normalized phasor profile in the complex plane. 312.4 Estimated cfin(ω) for phantom No.3. Plotted in dark solidline is the selected region for estimation of average speed. Atlower frequencies no peak occurs inside the block length. . . . 322.5 (a) Estimated value of the Young’s modulus and (b) relax-ationtimesforthefivephantomsatdifferentfrequenciesofex-citation. The lower frequency graphs (4–40Hz) are from con-ventional rheometry and the higher frequency graphs (above40Hz) are from model fitting. . . . . . . . . . . . . . . . . . . 342.6 (a) The result of the excitation with a plate which covers thetop surface of the phantom (b) The result of the excitationwith a point source. . . . . . . . . . . . . . . . . . . . . . . . 362.7 Obtained Young’s modulus as a function of the volume por-tion of Plastic Softener or Hardener added during the fabri-cation of the phantoms . . . . . . . . . . . . . . . . . . . . . . 37viiiList of Figures3.1 Infinite medium; shear beam patterns for three different fre-quencies 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 100Hzrespectively. Note that only a cubic portion of the mediumnear the axis of symmetry is shown; The medium extends inall directions to infinity. . . . . . . . . . . . . . . . . . . . . . 623.2 Infinite medium; shear beam patterns for three different fre-quencies 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 100Hzrespectively. Note that only a cubic portion of the mediumnear the axis of symmetry is shown; The medium extends inall directions to infinity. . . . . . . . . . . . . . . . . . . . . . 633.3 The first four modes of the cylindrical bar . . . . . . . . . . . 643.4 Permissible beam patterns (modes) in an infinite cylindricalrod for waves having a frequency of 100Hz. Three (and onlythree) phase speeds are possible. (a), (b), (c) radial com-ponents of the displacement (d), (e), (f) z-component of thedisplacement for phase speeds of 3.0339m/s, 3.719m/s, and5.630m/s respectively. . . . . . . . . . . . . . . . . . . . . . . 664.1 Displacement vector in the cylindrical coordinate system. . . 854.2 The magnitude of the displacement in a 1D model of thewave equation with the purely elastic wave speed of 3.16m/s,length of 6cm, and Young’s modulus of 10kPa. (a) and (b)are the magnitude and contour plots for the case in whichthe viscosity is 5Pas. The frequencies denoted by case (i)and case (ii) are almost equal. (c) and (d) show the the plotsfor the case that the viscosity of the medium is 30Pas. Inthis case, f1 over-estimates the actual value, while f2 resultsin a relatively small under-estimation. . . . . . . . . . . . . . 904.3 The geometry of the 3D region used for FEM simulations. . . 924.4 The magnitude (a) and the phase (b) of the transfer func-tion between the center point of the phantom and the middlepoint at the top. The solid line shows the case where thefour vertical sides of the phantom are free to move and thedashed lines are for the case where the sides of the phantomare confined from bulging. The gray dash-dot line shows thetransfer function for a one-node excitation at the middle point. 94ixList of Figures4.5 Effect of the viscosity on the magnitude of the transfer func-tion between the center point and the middle point. The ex-citation was applied to the entire top surface while the sideswere free to bulge. . . . . . . . . . . . . . . . . . . . . . . . . 954.6 Transient step response for different boundary conditions . . 1024.7 Wave propagation in a soft PVC phantom with E = 17kPaand an axial length of 4cm under two different boundary con-ditions when a step excitation is applied. (a) and (c) showthe displacements of tissue at different depths for a flat anda point-source excitation respectively. The location of thewavefront peaks are shown versus time for these two bound-ary conditions in (b) and (d), where a line has been fitted tothe data to estimate the speed. The resonance frequency ofthe phantom is also a function of the elastic modulus, geom-etry and boundary conditions. . . . . . . . . . . . . . . . . . . 1035.1 Acquisition time-line for two frames of B-mode data on atwelve element probe. . . . . . . . . . . . . . . . . . . . . . . 1105.2 Acquisition time-line for two frames of high frame rate acqui-sition data on a twelve element probe . . . . . . . . . . . . . 1115.3 A hypothetical acquisition time-line for delay-free simultane-ous sampling on a twelve element probe . . . . . . . . . . . . 1135.4 Flow chart of sequence programming for high frame rate im-age acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . 1175.5 Acquisitiontime-lineinvolvingreprogrammingdelayonatwelveelement probe . . . . . . . . . . . . . . . . . . . . . . . . . . . 1185.6 The extension of a limited duration signal to a periodic signal 1205.7 Signal processing pipeline of the high frame rate ultrasoundsystem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1205.8 Schematics of (a) the side view and (b) the front view of theexperimental setup developed for testing the system. Notethe choice of the coordinate axes. . . . . . . . . . . . . . . . . 1235.9 Different steps in the phase correction process: the phase(a) before any compensation (b) after compensation for thedelays on each scan-line (c) after compensation for the delaysbetween the scan-lines in each sector (d) after compensationfor the delays between the sectors. . . . . . . . . . . . . . . . 125xList of Figures5.10 (a), (b), (c), (d) the amplitude, and (e), (f), (g), (h) corre-sponding phase of the displacement phasors at different fre-quencies of excitation. From left to right the frequencies are75, 100, 150, and 200Hz. . . . . . . . . . . . . . . . . . . . . . 1265.11 Spatial derivatives of the axial displacement phasor U(x,y) atthe excitation frequency of 100Hz. (a) |∂U/∂x| (b) |∂2U/∂x2|(c) |∂U/∂y| (d) |∂2U/∂y2| (e) ∠∂U/∂x (f) ∠∂2U/∂x2 (g)∠∂U/∂y (h) ∠∂2U/∂y2 . . . . . . . . . . . . . . . . . . . . . . 1275.12 Outline of the geometry of the inclusion phantoms used inthe elastacity experiments. . . . . . . . . . . . . . . . . . . . . 1285.13 (a) B-mode and (b) elastogram images of a phantom with acircular hard inclusion. . . . . . . . . . . . . . . . . . . . . . . 1295.14 (a) B-mode and (b) elastogram images of a phantom with acircular soft inclusion. . . . . . . . . . . . . . . . . . . . . . . 1296.1 Obtained Young’s modulus as a function of the volume por-tion of Plastic Softener or Hardener added during the fabri-cation of the phantoms . . . . . . . . . . . . . . . . . . . . . . 137A.1 Simple schematic of the rheometry device in a conventionalrheometry configuration. In this configuration the force sen-sor measures the force on the moving plate and moves withit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144A.2 Simple schematic of the rheometry device in a conventionalrheometry configuration. In this configuration the force sen-sor measures the force on the static plate and does not move. 145A.3 The infinitesimal volume which is used in the definition ofElastic modulus. . . . . . . . . . . . . . . . . . . . . . . . . . 146xiList of Abbreviations1D One Dimensional2D Two Dimensional3D Three DimensionalA-line Amplitude line, the amplitude of the Hilbert transform of an RF-lineARF Acoustic Radiation ForceARFI Acoustic Radiation Force ImagingB-mode Brightness mode, the most common type of medical ultrasound imagesDC motor Direct Current motorEC European CommissionECG EchoCardioGramEDTA EthyleneDiamineTetraAcetic acidFEM Finite Element Modeling or MethodMRE Magnetic Resonance ElastographyMRI Magnetic Resonance ImagingMSD Mass Spring Damper modelPC Personal ComputerPDE Partial Differential EquationPRF Pulse Repetition FrequencyPVA PolyVinyl AlcoholPVA-C PoylyVinyl Alcohol CryogelPVC PolyVinyl ChlorideP-wave Primary wave, defined in seismologyRF Radio Frequency, RF-data or RF-line is the signal received by an ultrasoundtransducerROI Region Of InterestS-wave Secondary wave, defined in seismologySNR Signal to Noise RatioTMM Tissue Mimicking MaterialxiiAcknowledgmentsThis work would not have been possible without the consistent supportof my parents, whose caring love has always been a source of motivationfor 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 liketo thank the members of the Robotics and Control Lab of University ofBritish Columbia Julian Guerrero, Orcun Goksel, Ehsan Dehghan, HaniEskandari, Reza Zahiri Azar, Sara Khosravi, and Jing Xiang for providing acollaborative and encouraging research environment. Living far away frommy family, my close friends Alireza Forghani, Soraya Kasnavi, and TriciaPang, have been the greatest source of support for me. In the end I shouldthank my trusted best friends Thomas Diego Prananta, Denis Tran, AndreLei, Budi Kartono, and Joanna Lam who have supported me with theirinvaluable friendship throughout these years. I wish all of them the best intheir lives and careers.xiiiStatement of Co-AuthorshipThis thesis is based on several manuscripts, resulting from collaborationbetween multiple researchers.A version of Chapter 2 appeared in the IEEE Transactions on Ultra-sonics, Ferroelectrics and Frequency Control in July, 2009. This article wasco-authored with Dr. Hani Eskandari, Dr. Tim Salcudean, and Dr. RobertRohling. The author’s contribution in that article was developing the mainidea, numerical simulation, phantom fabrication, experimental evaluationof the algorithms and writing the manuscript. Dr. Eskandari contributedto 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. Sal-cudean helped with their numerous suggestions in the course of the dataacquisition and analysis. They also assisted by editing the manuscript.A version of Chapter 3 appeared in the Journal of Acoustical Society ofAmerica in September, 2009. This article was co-authored with Dr. TimSalcudean and, Dr. Robert Rohling. The author’s contribution in the ar-ticle was developing the idea, mathematical analysis and proof, numericalsimulation, and writing the manuscript. Dr. Salcudean and Dr. Rohlingassisted by editing the manuscript.A version of Chapter 4 appeared in the “Physics in Medicine and Biol-ogy” in July, 2009. This article was co-authored with Dr. Hani Eskandari,Dr. Tim Salcudean, and Dr. Robert Rohling. The original idea that thelongitudinal waves can travel at a low speed was developed by the author.Dr. Eskandari made several contributions in particular adding viscosity tothe analysis, analyzing the solution for a point source placed on a cylin-der, and devising one of the methods for determining the wave speed fromdisplacement spectra. The second method was the idea of the author. Dr.Eskandari also performed the FEM simulations and the analysis of the re-sults. The experiments were first performed by the author and suggested toDr. Eskandari. However, the actual data used in the article were collectedby Dr. Eskandari. The manuscript was written by Dr. Eskandari. The au-thor was involved at different stages of the manuscript through discussionsxivStatement of Co-Authorshipwith Dr. Eskandari. Dr. Salcudean and Dr. Rohling assisted with theirsuggestions and editing the manuscript.A version of Chapter 5 has been submitted for publication. This arti-cle was co-authored with Alex Brant, Dr. Tim Salcudean and Dr. RobertRohling. The author’s contribution in that article was developing the mainidea, and designing and carrying out the experiments. Alex Brant convertedthe author’s Matlab code into C++ and incorporated it into a larger pro-gram running on the ultrasound machine for real-time data acquisition. Dr.Salcudean and Dr. Rohling assisted with their suggestions throughout thecourse of the research and editing the manuscript.xvChapter 1Introduction1.1 MotivationEach of the medical imaging modalities provides the clinicians with sometype of information about the internal characteristics of tissue. The infor-mation 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 scatter-ing properties of the tissue in the X-ray band of the electromagnetic wavespectrum. The contrast in magnetic resonance (MRI) images depends onthe characteristics of the hydrogen atom in the molecules in tissue, or itsconcentration. The contrast in ultrasound images is a function of the ab-sorption and scattering properties of the tissue for mechanical waves in the1-10MHz range of frequencies.Elasticity imaging or elastography refers to medical imaging techniqueswhich provide the clinicians with images and measurements of the viscoelas-tic properties of tissue [1–8]. Changes in the viscoelastic properties of tissueare often associated with pathology. In the diagnosis of breast tumors, thephysician often examines the breast by palpation in order to feel the pres-ence of hard lumps. The same is done for liver tumors or through digitalrectal examination for prostate cancer. As useful as it is, palpation is not anaccurate way of determining the tissue stiffness. Images of tissue elasticity,however, provide quantitative measures for the physician to make furtherdecisions for the diagnosis or treatment.1.2 BackgroundIn spite of their differences, all the elastography techniques developed so farshare a general methodology for obtaining tissue mechanical properties. Inthis methodology, the tissue is subjected to some form of deformation [9–17],as it is being imaged by a medical imaging system. The images obtainedfrom the deforming tissue are further processed to measure tissue motion,namely displacements or velocities [11,18–27]. The tissue motion is then1Chapter 1. Introductionused in an inverse problem to find its mechanical properties [11,28–40]. Theimages of tissue elasticity are also called elastograms in the literature.When the mechanical properties of a sample are known, and the sampleis subjected to forces, the motion of the material can be found by solvinga forward problem. Therefore the reverse of this process, that is of findingthe mechanical properties of the material from its displacements is calledsolving the inverse problem.The elastography techniques can be categorized based on the excita-tion technique, the medical imaging technique, and the inverse problemapproach.Type of ExcitationThe excitation can be static, or dynamic. A dynamic excitation can besingle frequency, or wideband. The excitation can be external, for instancea 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 ModalityMost often, the medical imaging technique used for elastography is eitherultrasound or magnetic resonance imaging (MRI). Each of these techniqueshas its own advantages and disadvantages. MRI can provide measurementsof tissue displacement with sub-micron accuracy and a high signal to noiseratio (SNR) [11,36,37,41,42]. It can measure displacements in all threespatial directions (x, y, and z), hence provides 3D vector displacements.Moreover the 3D vector displacements can be measured at all spatial loca-tions of a 3D volume. This is the most complete set of information aboutthe tissue displacement, as it completely describes the way the tissue is de-forming. However obtaining such MRI images can take from a few minutesto a portion of an hour. The motion of the patient during the prolongedacquisition interval reduces the SNR and the accuracy of the measurements.Moreover MRI is not as readily available as ultrasound and is more expensiveto perform.Ultrasound imaging devices on the other hand are relatively inexpensiveand more readily available [18–27]. The acquisition process is much faster,and with the improvements in the elastography techniques, real-time ultra-sound elastography is not a far-fetched goal anymore. The disadvantages arethat the measurements are less accurate than MRI and suffer from a lowerSNR. The most important shortcoming of ultrasound, however, is that in its2Chapter 1. Introductionmost basic form, it only measures 1D displacements in the direction of theultrasound wave propagation, and the measurements are carried out onlyover a 2D plane, namely the plane of imaging. Even if 3D measurementscould be made available through 2D or mechanical 3D arrays, ultrasoundresolution in the lateral direction has fundamental limits which are typicallyan order of magnitude less than the axial resolution. In this case, a com-plete description of the tissue motion is not available and the mechanicalproperties need to be estimated based on the partial information.Inverse Problem TechniquesThe inverse problem techniques used in elastography are quite versatile andto some extent, depend on the excitation technique used. One of the simplesttechniques used in manual ultrasound elastography is to show the inverseof the strain. In this technique the physician applies pressure on the tissueusing the ultrasound probe itself, and it is assumed that this pressure causesa uniform constant stress in the material. Under these conditions the strainbecomes inversely proportional to the tissue elasticity. Therefore in thiscase, 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 distributionof elasticity values is assumed for the nodes of the mesh. This could be a uni-form elasticity distribution for instance. The forward problem is then solvedbased on these elasticity values, using FEM techniques, and the expecteddisplacements are found. The expected displacements are then comparedto the measured displacements. As these would not match, an optimizationscheme 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 mechan-ical wave equation [11,13–16,23,36–38,46–50]. When a dynamic excitationis used, the motion of the tissue is governed by a wave equation. The formof the equation depends on the level of modeling used for the tissue, forinstance whether anisotropy or viscosity are taken into account or not. Foran isotropic purely elastic homogeneous medium, in the absence of body3Chapter 1. Introductionforces, the equations take the form:ρ∂2u∂t2 = (λ+µ)∂∆∂x +µ∇2u (1.1)ρ∂2v∂t2 = (λ+µ)∂∆∂y +µ∇2v (1.2)ρ∂2w∂t2 = (λ+µ)∂∆∂z +µ∇2w (1.3)The wave equation is a partial differential equation; the independent vari-ables are the spatial coordinates (x,y,z), and the time t. The dependentvariables are the 3D displacement components (u,v,w). The mechanicalproperties 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 andis defined in terms of the displacements by,∆ = ∂u∂x + ∂v∂y + ∂w∂z (1.4)From a control systems perspective, the inverse problem can be viewedas an identification problem; the measured quantities are the dependentvariables 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 ispossible to carry out a full inversion of the wave equation [41,42]. However,as is the case in ultrasound elastography, the 3D vector displacement isnot readily available. Moreover, in some magnetic resonance elastography(MRE) techniques, to decrease the imaging time, not all three componentsof the displacements are measured. In these cases, the inversion techniquesuse 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 thephase speed, then estimate the tissue elasticity [13–16,36–38,47,48,50]. Thephase speed is the speed at which equiphase surfaces propagate throughthe medium. Under certain excitation and boundary conditions, the phasespeed cphase for shear waves is related to the shear modulus µ,cphase =radicalbiggµρ . (1.5)4Chapter 1. IntroductionDifferent techniques exist for estimating the phase speed. In transientelastography, the excitation is in the form of a short pulse [13–15]. Thepropagationspeedofthispulseisthenmeasured. Insome ofthese techniquesthe measured quantity is exactly the phase speed.In another class of dynamic elastography techniques, the excitation iscontinuously applied to the tissue, and the measured displacements are inthe steady state [11,36,37,48]. In these techniques the displacement areanalyzed in the phasor domain. Various techniques such as algebraic inver-sion and phase gradient have been proposed for estimating the phase speed.The underlying problem which all these techniques strive to solve is that ofestimating 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 availableover a 3D volume, a complete inversion can be carried out. However in thecase of partial displacement data, where the phase speed is estimated, a rela-tionship exists between the resolution of the elastograms and the frequencyof excitation. This is due to the fact that to measure the local frequency ofa signal with a higher accuracy, more of its periods should be observed. Asa result, the higher the frequency of excitation, the shorter the wavelengthof the waves, and the higher the resolution of the elastograms. Elastogramswith higher resolution are clinically preferred, because they can show smallertumors and calcifications and lead to earlier diagnosis. But the higher fre-quencies are attenuated more in real tissue, so there is a practical limit tothe number of wavelengths that can be observed in a given ROI. Neverthe-less, there are still limitations on the ability of the ultrasound machine tocapture motion at these frequencies.1.3 Thesis ObjectivesThis thesis studies the inversion problem of elastography under the followingconditions:• 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 waveequation-based approach can be used in ultrasound elastography to pro-vide images of tissue elasticity. To examine this hypothesis, the following5Chapter 1. Introductionobjectives are set:• Finding a method for accurately measuring the elasticity of tissue andtissue mimicking material samples for setting gold standards.• Fabricating tissue mimicking materials with controllable elasticity.• Establishing the relationship between the phase speed and tissue elas-ticity through a theoretical analysis of the elastic wave equation.• Developing an ultrasound elastography system which uses the devel-oped algorithm for forming elastograms.The following contributions have been made towards achieving the ob-jectives:• A novel rheometry method has been devised which is capable of ac-curate measurements of tissue elasticity and viscosity on test blocks.The accuracy of the device has been shown by a comparison of theresults with conventional rheometry.• Poly Vinyl Chloride (PVC) tissue mimicking materials have been fab-ricated and their elasticity characterized based on the volume portionof added softener or hardener.• The limitations of the shear wave speed equation (1.5) have been es-tablished. In particular experimental data have been provided forwhich equation (1.5) does not apply.• It has been established theoretically that a universal relationship doesnot exist between the phase speed and tissue elasticity. However underspecial boundary and excitations, a relationship can exist.• The relationship between the phase speed and tissue elasticity is es-tablished for small and large profile exciters placed on finite mediaunder different boundary conditions, by using FEM simulations andalso experimentally.• A novel high-frame-rate ultrasound system has been devised whichuses the idea of sector subdivision to achieve the increased frame rate.Unlike other systems that use sector subdivision, the devised systemcompensates for the time delays introduced in the acquisition processby a novel frequency domain phase compensation scheme.• Conditions have been found which should be met by the mechanicalexcitation in order for the high-frame-rate system to work properly.• The high-frame-rate system has been used in conjunction with a phasespeed based inversion algorithm to develop and ultrasound elastogra-phy system. The performance of the elastography system has beentested on tissue mimicking material.6Chapter 1. Introduction1.4 Chapter SummaryThe overall goal of this thesis is to study the application of wave equation-based inversion techniques in ultrasound elastography. This thesis has beenwritten in the manuscript-based style, as permitted by the Faculty of Grad-uate Studies at the University of British Columbia. In the manuscript-basedthesis, each chapter represents an individual work that has been published,submitted or prepared for submission to a peer reviewed publication. Eachchapter is self-sustained in the sense that it includes an introduction to thework presented in that chapter, the methodology, results and discussion.The references are summarized in the bibliography found at the end of eachchapter. The appendices pertaining to each chapter are presented at theend of the thesis.In Chapter 2 we show that for a finite block of material, the phase speedfor longitudinal waves could be related to the shear modulus as in,cphase =radicalbigg3µρ . (1.6)Furthermore, based on this relationship, a novel method is introduced forrheometry, i.e. determining the elasticity and viscosity of sample blocks ofmaterial. This method is more suitable for studying soft tissue specimensthan conventional rheometry. The accuracy of the method is tested againstconventional rheometry.In Chapter 3 a further theoretical analysis is carried out on the phasespeed which reveals that neither (1.5) nor (1.6) is a universal relation. Infact, for the same material, with the same shear modulus, the phase speedcan take different values, depending on the beam form of the wave, and theboundary conditions.In Chapter 4 the results of Chapter 3 are corroborated by using FEMand experiments. With a 3D FEM simulation, it is shown that a pointexciter creates longitudinal waves which travel at the phase speed given by(1.5) but a large profile exciter creates longitudinal waves which travel atthe phase speed given by (1.6). Both harmonic and transient analyses areperformed. Similar results are obtained from phantom experiments.In Chapter 5 a method is presented for acquiring displacement phasorsat a high-frame-rate. The low frame rate of conventional ultrasound systemslimits the scope of motions they can track. Namely, the Nyquist criterionimposes a limit on the highest frequency of excitation that can be used withthese systems. This limit is in the range of 50Hz to 200Hz, depending on7Chapter 1. Introductionthe depth of imaging and field of view. In the current setting, to form elas-tograms of higher resolution it is necessary to excite the tissue at higherfrequencies, which in turn requires a higher frame rate on the imaging sys-tem. The proposed method can easily be implemented on most commercialultrasound systems without the need for a hardware upgrade. The effective-ness of the system in forming elastograms is demonstrated experimentallyon tissue mimicking material with soft and hard inclusions.In Chapter 6 the thesis is concluded by a brief summary and a discussionwhich shows the place of the contributions of this thesis in the bigger pic-ture of tissue elastography. Recommendations are made for future researchdirections.Two appendices are appended to the thesis. In Appendix A the limita-tions of conventional rheometry are analyzed. In Appendix B a literaturesurvey of the most commonly used tissue mimicking materials is presentedtogether with a section on fabricating Poly Vinyl Chloride phantoms. InAppendix C the method for deriving the phasor of the displacements is de-scribed together with two typical inversion techniques based on the waveequation.8References[1] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastog-raphy: A quantitative method for imaging of elasticity of biologicaltissues,” Ultrasonic Imaging, vol. 13, pp. 111–134, 1991.[2] K. Parker, L. Gao, R. Lerner, and S. Levinson, “Techniques for elas-tic imaging: A review,” IEEE Engineering in Medicine and BiologyMagazine, vol. 15, no. 6, pp. 52–59, 1996.[3] L. Gao, K. J. Parker, R. M. Lerner, and S. Levinson, “Imaging of theelastic 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 andGases. Academic Press, 2001, vol. III, ch. Elastic Properties of SoftTissues.[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 ultra-sound,” Journal of Medical Ultrasonics, vol. 29, pp. 155–171, 2002.[7] J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imag-ing elastic properties of biological tissues,” Annual Review BiomedicalEngineering, vol. 5, pp. 57–78, April 2003.[8] K. J. Parker, L. S. Taylor, S. Gracewski, and D. J. Rubens, “A unifiedview 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 char-acterization and echographic imaging,” in Proceedings of the 7th Euro-pean Communities Workshop, 1987, nijmegen, Netherlands.9Chapter 1. Introduction[10] R. Lerner, K. Parker, J. Holen, R. Gramiak, and R. Waag, “Sono-elasticity: Medical elasticity images derived from ultrasound signals inmechanically vibrated targets,” Acoustical Imaging, vol. 16, pp. 317–327, 1988.[11] R. Sinkus, J. Lorenzen, D. S. amd M. Lorenzen, M. Dargatz, andD. Holz, “High-resolution tensor mr elastography for breast tumor de-tection,” Physics Med. Biol., vol. 45, pp. 1649–1664, 2000.[12] R. Souchon, L. Soualmi, M. Bertrand, J.-Y. Chapelon, F. Kallel, andJ. Ophir, “Ultrasonic elastography using sector scan imaging and aradial compression,” Ultrasonics, vol. 40, pp. 867–871, 2002.[13] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulusimaging with 2-d transient elastography,” IEEE Trans. Ultrason. Fer-roelect. 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 elas-tography in anisotropic medium: Application to the measurement ofslow 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 imag-ing: 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 andmotion 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-adaptivemotion tracking of ultrasound image sequences using a deformablemesh,” IEEE Trans. Med. Imaging, vol. 17, no. 6, pp. 945–956, 1998.[19] E. Konofagou, T. Varghese, J. Ophir, and S. Alam, “Power spectralstrain estimators in elastography,” Ultrasound Med. Biol., vol. 25, no. 7,pp. 1115–1129, 1999.10Chapter 1. Introduction[20] E. E. Konofagou and J. Ophir, “Precision estimation and imaging ofnormal 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, “Di-rect 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 IEEEUltrasonics Symposium, 2002, pp. 1811–1820.[23] J. Bercoff, R. Sinkus, M. Tanter, and M. FinkJ, “3d ultrasound-baseddynamic and transient elastography: First in vitro results,” in IEEEUltrasonics Symposium, 2004, pp. 28–31.[24] K. Hoyt, F. Forsberg, and J. Ophir, “Investigation of parametric spec-tral estimation techniques for elasticity imaging,” Ultrasound Med.Biol., vol. 31, no. 8, pp. 1109–1121, 2005.[25] ——, “Comparison of shift estimation strategies in spectral elastogra-phy,” Ultrasonics, vol. 44, pp. 99–108, 2006.[26] ——, “Analysis of a hybrid spectral strain estimation technique in elas-tography,” Phys. Med. Biol. 51, vol. 51, pp. 197–209, 2006.[27] J.-L. Gennisson, T. Deffieux, R. Sinkus, P. Annic, M. Pernot, F. Cud-eiro, G. Montaldo, M. Tanter, M. Fink, and J. Bercoff, “A 3d elastogra-phy system based on the concept of ultrasound-computed tomographyfor 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 noninvasivedetermination of material parameters from a knowledge of elastic dis-placements: 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 modu-lus 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 elas-ticity imaging of soft tissues,” IEEE Transactions on Nuclear Science,vol. 41, no. 4, pp. 1639–1648, August 1994.11Chapter 1. Introduction[31] F. Kallel and M. Bertrand, “Tissue elasticity reconstruction using linearperturbation 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 re-construction 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 inverseproblems in elasticity imaging using the adjoint method,” Inverse Prob-lems, vol. 19, pp. 297–313, 2003.[35] D. Fu, S. Levinson, S. Gracewski, and K. Parker, “Non-invasive quan-titative reconstruction of tissue elasticity using an iterative forward ap-proach,” Phys. Med. Biol. 45 (2000) 14951509, vol. 45, pp. 1495–1509,2000.[36] A. Manduca, T. Oliphant, M. Dresner, J. Mahowald, S. Kruse, E. Am-romin, J. Felmlee, J. Greenleaf, and R. Ehman, “Magnetic resonanceelastography: Non-invasive mapping of tissue elasticity,” Medical ImageAnalysis, vol. 5, pp. 237–2540, 2001.[37] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complex-valued stiffness reconstruction by for magnetic resonance elastographyby algebraic inversion of the differential equation,” Magnetic Resonancein Medicine, vol. 45, pp. 299–310, 2001.[38] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelka-ram, and J. Culiolic, “Measurement of viscoelastic properties of ho-mogeneous soft solid using transient elastography: An inverse problemapproach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec2004.[39] B. Robert, R. Sinkus, B. Larrat, M. Tanter, and M. Fink, “A newrheological 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 withsurface wave speed (l),” J. Acoust. Soc. Am., vol. 122, no. 5, pp. 2522–2525, November 2007.12Chapter 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 mrelastography,” Magnetic Resonance Imaging, vol. 23, pp. 159–165, 2005.[42] R. Sinkus, M. Tanter, S. Catheline, J. Lorenzen, C. Kuhl, E. Son-dermann, and M. Fink, “Imaging anisotropic and viscous propertiesof breast tissue by magnetic resonance-elastography,” Magnetic Reso-nance in Medicine, vol. 53, no. 2, pp. 372–387, 2005.[43] A. A. Oberai, N. H. Gokhale, M. M. Doyley, and J. C. Bamber, “Eval-uation 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 moduluselastography,” Ultrasound Med. Biol., vol. 31, no. 6, pp. 787–802, 2005.[45] M. M. Doyley, S. Srinivasan, E. Dimidenko, N. Soni, and J. Ophir, “En-hancing the performance of model-based elastography by incorporatingadditional a priori information in the modulus image reconstructionprocess,” 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 us-ing 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 andelastic properties ofsoft tissue using supersonic shear imaging,” in IEEEUltrasonics Symposium, 2003, pp. 925–928.[48] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity andviscosity 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. Mon-grain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel me-chanical properties with four ultrasound elastography methods andcomparison with gold standard testing,” IEEE Trans. Ultrason. Fer-roelect. Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007.[50] A. Romano, J. Bucaro, P. Abraham, and S. Dey, “Inversion methodsfor the detection and localization of inclusions in structures utilizing13Chapter 1. Introductiondynamic surface displacements,” in Proceedings of SPIE Vol. 5503, vol.5503, 2004, pp. 367–374.[51] E. Turgay, S. Salcudean, and R. Rohling, “Identifying mechanical prop-erties of tissue by ultrasound,” Ultrasound Med. Biol., vol. 32, no. 2,pp. 221–235, 2006.[52] H. Eskandari, S. Salcudean, and R. Rohling, “Viscoelastic parameterestimation based on spectral analysis,” IEEE Trans. Ultrason. Ferro-elect. Freq. Contr., vol. 55, no. 7, pp. 1611–1625, 2008.[53] B. Boashash, “Estimating and interpreting the instantaneous frequencyof 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 asignal – part 1:fundamentals,” Proceedings of IEEE, vol. 80, no. 4, pp.540–568, April 1992.14Chapter 2Measurement of ViscoelasticProperties of TissueMimicking Material UsingLongitudinal WaveExcitation∗2.1 IntroductionThe viscoelastic characterization of tissue and tissue mimicking material isof interest in many areas of biomedical research. The variability in tissuestiffness due to anatomical features, as well as the changes incurred in tissuestiffness due to pathology, has led to the development of elastography as aclinical imaging modality [1–5]. In general, elastography involves observa-tion of mechanical deformation patterns inside tissue by means of a clinicalimaging system and use of these patterns to infer mechanical properties ofthe 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 etal. proposed quasi-static compression [8]. Rudenko et al. used acoustic ra-diation force of focused ultrasound to create shear waves inside tissue [9,10].Sandrin et al. used external shakers which create longitudinally polarizedshear waves [11–13]. Excitation by two or more sources of slightly differentfrequencies has also been considered. The sources are focused ultrasound∗A version of this chapter has been published: A. Baghani, H. Eskandari, S.E. Salcud-ean, R. Rohling, “Measurement of Viscoelastic Properties of Tissue Mimicking MaterialUsing Longitudinal Wave Excitation,” IEEE Trans. on Ultras. Ferr. Freq. Cont., vol. 56,no. 7, pp. 1405-1418, Jul. 200915Chapter 2. Elastography Using Longitudinal Wavesin the work by Fatemi and Greenleaf [14]. They are shear wave sources inthe work by Hoyt et al. [15]. Among these works, only Ophir et al. usedcompressional excitation, and only in the quasi-static settings. In all theother cases, either dominantly shear excitation is used, or the compressionalexcitation is ignored or suppressed by taking the curl of the displacementfield.To quantitatively study the effectiveness of elastography imaging meth-ods in providing images of viscoelastic properties of tissue, it is of interest totest these methods on phantoms with known anatomical features of knownviscoelastic properties. A number of researchers are involved in the devel-opment of phantoms with adjustable viscoelastic properties [16–23].Some researchers have used conventional rheometry to validate the re-sults obtained from their elastography techniques [24,25]. Other researchershave developed specific methods for measuring viscoelastic properties of thephantom and soft tissue specimens [26,27]. Srinivasan et al. used a nanoin-denter to measure elasticity [28] . Kalanovic et al. compared the results ofelasticity measurements with a commercial indenter to the results obtainedfrom their developed shear rheometry device [29]. Egorov et al. have de-veloped a soft tissue elastometer which uses a small profile indenter and ascale to measure the applied force [30]. They concluded that the validity oftheir results depends on the size of the tissue sample relative to the size ofthe indenter. Samani et al. also used a small profile indenter to measureelasticity [31]. FEM was used in their work to study the geometric constantrelating the force to the indentation. In their work magnetic resonanceimaging based elastography was also proposed for the purpose of rheometry.Mechanical wave propagation in elastic continua has been studied ex-tensively and is well understood [32,33]. Some viscoelastic continua havealso been studied with certain assumptions about the viscoelastic behaviorof the matter [34,35]. In general the stress-strain relationship for a materialcan be very complex. The assumption usually made is that the stress-strainrelationship 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]. Longitudi-nally polarized transient shear waves were produced mechanically and theresulting displacements were tracked by ultrasound speckle tracking to givemeasurements of the wave speed and attenuation through an inverse prob-lem approach. It was shown that a Voigt model fits the visoelastic behaviorof the living tissue better than a Maxwell model.It is the objective of this research to:16Chapter 2. Elastography Using Longitudinal Waves• Show theoretically and experimentally that low frequency longitudinalwaves can propagate in a tissue mimicking material of finite dimen-sions at speeds which are much lower than those normally associatedwith ultrasound propagation.• Show experimentally that it is possible to use ultrasound-based speckletracking techniques to track the longitudinal waves.• Show theoretically and experimentally that the measured displace-ments caused by longitudinal wave propagation can be used to esti-mate the viscoelastic properties of the tissue mimicking material.To this end, a new framework is developed for viscoelastic characterizationof tissue and phantom specimens by using longitudinal wave excitation andvalidated against the conventional rheometry.We start by a brief review of longitudinal wave equations in elastic mediain Section 2.2. The experimental setup, material, and methods used arepresented in Section 2.3. Section 2.4 is devoted to the presentation of ourrheometry methods. The experimental results are presented in Section 2.5and discussed in Section 2.6. We conclude the paper in Section 2.7.2.2 Theory of Wave MotionIn this section the theory of wave motion in elastic media is briefly reviewedand the notation used in the article presented. The presentation is adoptedfrom [33]. We start with the model of an elastic medium. We use u(x)to denote the displacement vector field. τij and epsilon1ij are used to denote thestress and strain tensors. By definition:epsilon1ij = 12parenleftbigg ∂∂xj ui +∂∂xiujparenrightbigg. (2.1)The fundamental conservation equations in a medium can then be expressedas:mass : ∂∂tρ+summationtexti ∂∂xi(ρ˙ui) = 0, (2.2)momentum : summationtextj ∂∂xj τji +ρfi = ρ¨ui , (2.3)angularmomentum : τij = τji i negationslash= j, (2.4)energy : ρ˙e = summationtextij τji˙epsilon1ij . (2.5)where ρ(x) is the density, f(x) is the force per unit mass vector field, ande(x) is the internal energy per unit mass. These equations hold regardless of17Chapter 2. Elastography Using Longitudinal Wavesthe type of the medium. To model the specific behavior of different media,they are concatenated by a set of constitutive equations. These equationsmodel the stress-strain behavior of the medium.Ifweassumeahomogeneousisotropicelasticmedium, undersmallstrains,the constitutive equations widely used areτij = λ(epsilon111 +epsilon122 +epsilon133)δij +2µepsilon1ij , (2.6)where δij is equal to one if i = j, and zero otherwise. λ and µ are called theLam´e constants.A number of other constants are also defined based on the Lam´e con-stants,E = µ(3λ+2µ)/(λ+µ) = Youngprimesmodulus, (2.7)ν = 12λ/(λ+µ) = Poissonprimesratio, (2.8)K = λ+ 23µ = bulkmodulus. (2.9)The inverse of the equation (2.6) can be expressed in terms of these constantsasepsilon1ij = 1+νE τij − νE(τ11 +τ22 +τ33)δij . (2.10)At this point, it is worth defining the scope of our problem in terms ofthese parameters in dealing with soft human tissue and tissue mimickingmaterial. It has been found experimentally that these materials are nearlyincompressible [37,38]. Hence the Poisson’s ratio ν for these materials is veryclose to 0.5. In this respect, the material under study is closer to liquidssuch 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 tosolids than liquids. A number of consequences follow from these peculiarproperties:ν ≈ 0.5, (2.11)λ greatermuch µ, (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.18Chapter 2. Elastography Using Longitudinal WavesThe equations (2.2) to (2.6) form a complete set of equations for thedisplacement fields, stresses and density and can be used to derive the mostgeneral form of the wave equation [32,33]. For our purposes it suffices toconsider some special cases. We are particularly interested in longitudinal(compressional) waves. The two cases of interest are plane waves in aninfinite medium, and plane waves in a finite block of medium.2.2.1 Longitudinal Plane Waves in an Infinite MediumIn 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 theplane 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 lon-gitudinal plane wave propagating in the x1 direction in an infinite elasticmedium,(λ+2µ) ∂2∂x21u1 = ρ∂2∂t2u1 . (2.15)This equation predicts that the wave speed would becinf =radicalBiggλ+2µρ . (2.16)Because of the assumptions, in this case the strain tensor takes the formepsilon111 = ∂∂x1u1 , (2.17)epsilon1ij = 0 otherwise. (2.18)The stress tensor can be found using the constitutive equations (2.6),τ11 = (λ+2µ)epsilon111 , (2.19)τ22 = λepsilon111 , (2.20)τ33 = λepsilon111 , (2.21)τij = 0 otherwise. (2.22)It is evident from these equations that an infinitesimal cubic volume in themedium experiences normal forces acting on all of its surfaces as the planewave propagates through that small volume.Before we move on to the next case, it is important to mention that19Chapter 2. Elastography Using Longitudinal Wavesthe wave speed which is widely reported and accepted in the literature forultrasound iscUS =radicalBiggKρ =radicalBiggλ+ 23µρ (2.23)where K is the adiabatic compressibility in the case of fluids and bulk mod-ulus in the case of solids [39]. If the tissue is modeled as a fluid and thelongitudinal wave propagation as an adiabatic phenomenon, the wave speedwould be cUS. If it is modeled as an elastic solid, the wave speed would becinf of equation (2.16). However the difference between these two numbersis much smaller than the accuracy of the practical measurements in thisresearch, so cinf will be used throughout this paper as the ultrasound speed.Based on ultrasound theory and a large body of evidence, when a smallsource of longitudinal vibrations (piezoelectric crystal) with sub-millimeterdimensions is placed on the skin, with the human body as the propagationmedium, the pulses travel with a speed given by (2.23) or (2.16). Thereforefor a sub-millimeter wave source of high frequency, the body acts as aninfinite medium. As will be shown in the next subsection, the same is nottrue for all cases of longitudinal wave propagation in solids.2.2.2 Longitudinal Plane Waves in a Finite MediumAnother case of special interest is the propagation of longitudinal waves inmaterial blocks. In this case the effects of finiteness of media should alsobe taken into account. More specifically, if a rectangular block of materialis subjected to a spatially uniform, possibly time varying force on one of itssurfaces, a longitudinal wave would propagate in the block.The study of waves of this type could be facilitated by making someassumptions about the stress and strain patterns inside the block. Twomain assumptions are made here: that the “strain is uniform on each crosssection of the block at every instant”, and that “plane sections remain planeunder stress”. These assumptions make it possible to model the phenomenonas a one-dimensional problem.We further elaborate on these assumptions. A medium with a finite crosssection acts as a waveguide for elastic waves. At each frequency of excitationa number of spatial modes is excited. Each of these modes has its ownwave speed associated with it. As a result the wave speed depends on thefrequency of excitation. Now consider an infinitely long rod of finite crosssection which experiences a spatially uniform low frequency compressionover its cross section. Under these conditions the dominant mode of the rod20Chapter 2. Elastography Using Longitudinal Waveshas a wave pattern which satisfies the assumptions made here [32]. Moreprecisely, at low frequencies only the first mode is excited and the amplitudeof vibrations is almost constant over the cross section. This justifies the useof the assumptions made here and the use of a 1D model.Now assume that the wave is propagating in the x1 direction. The factthat the block has free boundaries on its sides causes the lateral stresses tovanish on the surfaces,τ22 = 0, (2.24)τ33 = 0, (2.25)τij = 0 i negationslash= j. (2.26)Because of the assumptions made, these conditions hold throughout theblock. If not, the infinitesimal volumes located around the boundary of theblock, would be (axially) compressed more, because of the absence of lateralstresses, than internally located infinitesimal volumes. As a result planesections wouldnotremainplane anymore whichcontradicts the assumptions.Moreover, the top and bottom boundaries of the block should be frictionlessslip boundaries to prevent any external lateral forces on them.Now the constitutive equation (2.10) givesepsilon111 = 1+νE τ11 − νEτ11= 1Eτ11 , (2.27)or equivalentlyτ11 = Eepsilon111= µ(3λ+2µ)λ+µ epsilon111 . (2.28)As a result, the equation of a longitudinal wave propagating in the x1direction in the block would beµ(3λ+2µ)λ+µ∂2∂x21u1 = ρ∂2∂t2u1 . (2.29)This equation predicts that longitudinal waves propagate in the block at a21Chapter 2. Elastography Using Longitudinal Wavesspeed ofcfin =radicalBiggµ(3λ+2µ)λ+µρ =radicalBiggEρ . (2.30)A comparison of (2.30) with (2.16) shows that the propagation speed oflongitudinal (compressional) waves is different in finite blocks as comparedto infinite media of the same material. For nearly incompressible materialcfin ≈ radicalbig3µ/ρ = √3cS where cS = radicalbigµ/ρ is the shear wave speed in theinfinite medium. This speed is in the order of 1 to 10m/s for soft tissuewhich is much smaller than cinf which is approximately 1540m/s.Since cfin is much smaller than cinf, it is possible to track the longitudinalwave propagation in blocks of tissue using pulsed ultrasound. Our objectiveis to perform rheometry measurements by studying the longitudinal wavepatterns inside finite blocks of tissue phantoms. In the next section theapparatus and methods used to create these wave patterns and collect dataare presented. The proposed rheometry methods for analyzing the collecteddata follow in Section 2.4.2.3 Experimental Methods and SetupIn this section the experimental setup is discussed, and the acquired data aredescribed. For ultrasound elastography studies, a device has been developedwhich is capable of applying dynamic compressions to a phantom while itis being imaged by an ultrasound probe. The detailed description of thisdevice can be found in [40].Figure 2.1 shows a simplified schematic of this device. The device hastwo plates between which the phantom is compressed. The upper plate isconnected to the chassis and does not move. This plate has an opening foran ultrasound probe, and fitted by a linear array ultrasound probe. Thelower plate is confined to move in one direction by a linear bearing. Thisplate is connected to a cam follower mechanism sitting on a cam driven bya controlled DC motor.For testing, a rectangular phantom block is placed between the two platesand the height of the upper plate is adjusted to have some pre-compressionon the phantom to ensure contact with the ultrasound transducer is main-tained. By (2.24) - (2.26) it is required that the top and bottom surfacesof the block only experience forces in the vertical direction. To match theexperimental conditions to the theoretical ones, ultrasound gel was used asa lubricant on the top and bottom surfaces of the phantom in contact with22Chapter 2. Elastography Using Longitudinal WavesFigure 2.1: Simple schematic of the rheometry devicethe 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 phan-tom material. The stiffness of this material can be adjusted by addingPlastic SoftenerTM or Plastic HardenerTM provided by the same company.SigmaCellr© Cellulose, Type 50 (Sigma-Aldrich Inc., St Louis, MO, USA)was added to the phantoms as the ultrasound scatterer to develop the specklepattern needed for ultrasound-based motion tracking. For each 100mL ofliquid plastic, 1mL of cellulose was added. We used a 40mm×30mm×30mmmold to fabricate five homogeneous blocks with different stiffnesses. The fab-rication parameters, together with the measured densities and ultrasoundspeeds cinf for these phantoms are summarized in Table 2.1. The densi-ties were measured by weighing the phantoms with a precision scale, andmeasuring the block dimensions with calipers and dividing the mass bythe calculated volume. The ultrasound speed was found by first measur-ing the depth to the bottom of the phantom on a B-mode ultrasound im-age which is displayed on the ultrasound machine assuming an ultrasoundspeed of 1540m/s, and then correcting the speed so that the measured depthmatches the height of the phantom. As can be seen, the density and cinfare not greatly affected by changing the proportion of Plastic to Softener orHardener.In all the experiments the lower plate was excited with Gaussian white-23Chapter 2. Elastography Using Longitudinal WavesTable 2.1: Fabrication and measured properties of the five test phantomsPlastic Hardener Softener density ultrasoundphantom (volume (volume (volume (kg/m3) speedNo portion) portion) portion) ±5% (m/s)±1%1 4 0 2 1020 14202 5 0 1 1030 14203 6 0 0 990 14104 5 1 0 1010 14405 4 2 0 1030 1410Table 2.2: Parameters of the ultrasound imaging systemparameter valueprobe type L9–4 linear arraycenter frequency 4MHzsampling frequency 20MHzdepth of focus 25mmimaging depth 40mmDoppler pulse repetition frequency 2.5kHznumber of samples acquired 75,472noise, band-limited to 4–150Hz, while RF-lines were collected for 30 secondsat 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 (UltrasonixMedical Corp., Richmond, BC, Canada). The parameters of the imagingsystem are summarized in Table 2.2.The acquired RF-line data were processed off-line to obtain the displace-ment profile of each point along the scan line. Since the models used are 1Dmodels, only the axial displacements were calculated. Lateral displacementsalso occur which could be used to improve the estimations in future work.The algorithm used for displacement calculations is described in [41]. Theparameters 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 timet. It is worth mentioning that both x and t are discrete variables in thisnotation.The discrete Fourier transform was then applied to um(x,t) to obtainthe phasor displacement profile Um(x,jω). Frequency domain averagingwas also used, as provided by the tfestimate function of MATLABr© to24Chapter 2. Elastography Using Longitudinal WavesTable 2.3: Parameters of the speckle tracking motion estimation algorithmparameter valuenumber of blocks 80block length 1mmblock overlap 50%displacement measurement relative (corresponding to block velocity)stretching not usedsmooth the data. The data were normalized at each frequency separatelyso that the element located at x = 40mm would have unit amplitude andzero phase. This normalization does not have any effect on the results, butit helps with the plotting of the phasor data.At each frequency ω, Um(x,jω) is a complex-valued function of x whichgives the steady-state displacement phasor of a point located at x inside thephantom, if the phantom were excited with frequency ω. A typical phasorprofile is plotted in Fig. 2.2. Parts (a) and (b) show the amplitude andphase of the phasor profile as a function of depth x, while part (c) shows thesame 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 viscoelasticproperties of the phantoms.To validate our estimation results, we collected conventional rheometrydata from the phantoms. This was achieved by a modification to the deviceof Fig. 2.1, where the ultrasound probe and upper plate were replaced bya force sensor (Load Cell MDB-2.5, Transducer Techniques, Temecula, CA,USA). The rheometry device was reported in a previous article [40]. Moredetails about the rheometry method used in this article can be found in theAppendix A.2.4 Proposed Rheometry MethodsIn this section two methods are introduced for analyzing the obtained datapresented in the previous section for estimation of viscoelastic properties oftissue phantoms. We start by deriving a viscoelastic model for the phantomblock to be used in these methods.Longitudinal wave propagation in a viscoelastic block can be modeled ina number of different ways [32,33]. The most complete model is obtainedby solving the three dimensional wave equation in a viscoelastic mediumwith boundary conditions adopted to match that of the block. Except for25Chapter 2. Elastography Using Longitudinal Waves(a)(b)(c)Figure 2.2: A typical normalized displacement phasor profile, belonging tophantom No.3 at a frequency of 90Hz. The (a) magnitude and (b) phaseof the phasor profile (c) the same normalized phasor profile in the complexplane. The profile is normalized to the phasor of the point at x = 40mm.26Chapter 2. Elastography Using Longitudinal Wavessome special cases for which analytical solutions have been derived [33], thismodel is quite difficult to solve and requires numerical methods such as finiteelement 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 firstcoordinate by x ≡ x1. In the frequency domain, this equation becomes(jω)2ρU = µ(3λ+2µ)λ+µ ∂2∂x2U= E ∂2∂x2U. (2.31)A next step in the modeling is the introduction of a first order viscosityterm. The Voigt model was shown to more closely follow the behavior ofthe soft tissue [36]. In this case the wave equation becomes(jω)2ρU = E ∂2∂x2U+Bjω ∂2∂x2U, (2.32)where B models first order viscosity effects. This is basically equivalent toassuming a Voigt model for the stress-strain relation of the material. Theratio B/E is called the relaxation time.Note that both of the equations (2.31) and (2.32) have the form(jω)2U = c2fin(ω) ∂2∂x2U , (2.33)where cfin(ω) is a complex frequency-dependent wave speed given by,c2fin(ω) =Eρ, for (2.31)Eρ +jBωρ , for (2.32) .(2.34)The steady state solution to (2.33) can be found easily. Consider thecharacteristic equation of this differential equation,s2 + ω2c2fin(ω) = 0, (2.35)27Chapter 2. Elastography Using Longitudinal Wavesand denote its roots by s1(ω) = −s2(ω) = α(ω) + jβ(ω). Then, the steadystate 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 resultU(0,jω) ≡ 0 = A1 +A2 , (2.37)and therefore A2 = −A1. The phasor solution (2.36) in this case could besimplified to obtain,U(x,jω) = 2A1 (cos(βx)sinh(αx)+jsin(β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 standingwave 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 proposetwomethodsfortheestimationofthe viscoelastic parameters ofthe phantomblocks.2.4.1 Peak to Node MethodThis method can be used to estimate Young’s modulus for a block of mate-rial, based on the longitudinal wave pattern observed inside the block. From(2.38) the amplitude of the solution, after some manipulation, can be foundas|U(x,jω)|2 = 4A21parenleftbigsinh2(αx)+sin2(βx)parenrightbig . (2.39)Now if the viscosity is small, the model (2.31) can be used and thesolution to the characteristic equation (2.35) would beα+jβ = j ωradicalBigEρ, (2.40)28Chapter 2. Elastography Using Longitudinal Wavesor,α = 0, (2.41)β = ωradicalBigEρ. (2.42)In this case, (2.39) becomes|U(x,jω)| = 2A1vextendsinglevextendsinglevextendsinglevextendsinglevextendsinglevextendsinglesin ωradicalBigEρxvextendsinglevextendsinglevextendsinglevextendsinglevextendsinglevextendsingle . (2.43)The distance of a node to a peak, dnp, of the wave pattern is equal to aquarter of the wavelength, thereforecfin = (4dnp) ω2pi , (2.44)E = ρc2fin . (2.45)In our case the node is at x =0mm and at each frequency, the peak is foundby searching for a local maximum in a filtered version of the amplitude ofUm(x,jω).In summary, the peak to node method estimates Young’s modulus basedon the peak to node distance of a measured phasor profile Um(x,jω) usingequation (2.45). But this method does not provide estimates of viscosity orrelaxation time. In the next subsection another method is presented whichalso estimates the viscosity.2.4.2 Model Fitting MethodA more general method to characterize the viscoelastic properties of thesample block is to use the characteristic equation (2.35) to determine thecomplex wave speed, cfin(ω), for different excitation frequencies, ω, experi-mentally, and then use a parameterization of the viscoelastic properties, suchas that presented in (2.34) to find the particular parameters of interest.Alternatively, one of the parameterized models (2.34) can be directlyfitted to the data. We present this approach here. The viscoelastic modelwas chosen and the general solution (2.36) was fitted to the data by leastsquares optimization over the parameters A1, A2, E and B. Thus for each29Chapter 2. Elastography Using Longitudinal Wavesfrequency the following cost function was minimized:J(ω) =summationdisplayx|Um(x,jω)−Um(x,jω;E,B)|2 . (2.46)Figure 2.3 shows a typical fitted model. As can be seen the model fitsthe 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 ResultsThe proposed rheometry methods of Section 2.4 were applied to the datacollected using the apparatus and methods of Section 2.3. The results forthe five phantoms follow in this section.2.5.1 Peak to Node ResultsThe peak to node method was used to analyze the measured phasor profilesUm(x,jω) for the five phantom blocks. Note that (2.44) and (2.45) could becalculated at any frequency for which there is a node and a peak present inthe phasor profile (see Fig. 2.2 (a)). As a result, this method cannot be usedat very low frequencies, since no peaks occur within the 40mm height of thephantom. Figure 2.4 shows the estimated cfin(ω) for the test phantom No. 3over the range of frequencies where this method could be used. The averagespeed over this frequency range, cfin, is equal to 7.5m/s and the standarddeviation is equal to 0.29m/s. The Young’s modulus calculated based onthis average speed is 55kPa.The same method was used to estimate longitudinal wave speed andYoung’s modulus for the other four test phantoms. The results are summa-rized in Table 2.4. Also presented are the values obtained from conventionalrheometry on the same blocks. As can be seen the results are generally ingood agreement and this supports the validity of the model and thus theassumptions.2.5.2 Model Fitting ResultsThemodelfittingapproachwasusedtofindthevaluesofYoung’sModulusEand relaxation time for the five phantom blocks. Figure 2.5 compares the re-sults to the conventional rheometry results. In Fig. 2.5(a) the model-fittingand rheometry values of the Young’s modulus, E, for different frequencies30Chapter 2. Elastography Using Longitudinal Waves(a)(b)(c)Figure 2.3: Fitted model to the measured displacement data. The data isthe 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 normalizedphasor profile in the complex plane.31Chapter 2. Elastography Using Longitudinal WavesFigure 2.4: Estimated cfin(ω) for phantom No.3. Plotted in dark solid lineis the selected region for estimation of average speed. At lower frequenciesno peak occurs inside the block length.Table 2.4: Estimated values of longitudinal wave speed cfin and Young’smodulus for the test phantoms using the peak to node methodpeak toaverage node conventionalwave freq. Young’s rheometryspeed range modulus Young’sphantom ±std used ±std modulus±10%No (m/s) (Hz) (kPa) (kPa) error1 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%32Chapter 2. Elastography Using Longitudinal WavesTable 2.5: Estimated values of longitudinal wave speed cfin and Young’smodulus for the test phantoms using the model fitting methodmodelfitting conventionalwave freq. Young’s rheometryspeed range modulus Young’sphantom ±std used ±std modulus±10%No (m/s) (Hz) (kPa) (kPa) error1 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 thenormalized mean squared error defined byNMSE =radicalBiggJ(ω)summationtextx|Um(x,jω)−mean(Um(x,jω))|2, (2.47)is better than 10%. This threshold was used to ensure that the model fittedthe data well enough for further conclusions to be valid. The rheometryplots are limited to 40Hz. As expected, the values of E from model fittingmatches with rheometry values. The average values and standard deviationsare summarized in Table 2.5.In Fig. 2.5(b) the optimized values of the relaxation times together withthe values obtained from the conventional rheometry are presented. Notethat both axes are in logarithmic scale in this figure. Again the data of theconventional rheometry are limited to 40Hz. As can be seen from this figure,the general trends and values of relaxation times from model fitting matchthat of conventional rheometry. Also a power-law decrease in the relaxationtime with the increase of the frequency is observed which is consistent withthe previous reports [40].2.6 DiscussionIt has been shown experimentally by Catheline et al. [36] that when a finiteblock of viscoelastic material is vibrated by a compressional exciter with a33Chapter 2. Elastography Using Longitudinal Waves(a)(b)Figure 2.5: (a) Estimated value of the Young’s modulus and (b) relaxationtimes for the five phantoms at different frequencies of excitation. The lowerfrequency graphs (4–40Hz) are from conventional rheometry and the higherfrequency graphs (above 40Hz) are from model fitting.34Chapter 2. Elastography Using Longitudinal Wavessmall cross sectional area, two waves propagate in the block: a compressionalwave which parts almost instantaneously at the speed of ultrasound, i.e.1540m/s, and a longitudinally polarized slow shear wave, which follows withthe shear wave speed, i.e. radicalbigµ/ρ which is around 1-10m/s. The shear wavecan be tracked by ultrasound speckle tracking and the displacements can beused to measure the viscoelastic properties of the block.We have shown here that when a finite block of viscoelastic material isvibrated by a compressional exciter with a large cross sectional area whichcovers the entire top surface of the block, instead of the longitudinally po-larized slow shear wave, a slow longitudinal wave propagates with the speedcfin = radicalbigE/ρ ≈ radicalbig3µ/ρ which is around 5-15m/s. This longitudinal wavecan also be tracked by ultrasound speckle tracking and the displacementscan be used to measure the viscoelastic properties of the block.To show how the size of the exciter affects the wave speed we conductedtwo different experiments with the same phantom block, namely block No.1, similar to the experiments of [36]. In the first experiment, a plate exciterwas used which covered the entire top surface of the phantom. The surfacewas lubricated with ultrasound gel to minimize the friction between the plateand the phantom. A single period of 50Hz sinusoid was used as a transientexcitation. The exciter was an open loop voice coil. Therefore the actualmechanical vibration transferred to the phantom was not exactly a singleperiod of 50Hz, but had some tail. One scan line (A-line) was sampled at5kHz and the displacements tracked. The displacements at different depthsare plotted as a function of time in Fig. 2.6(a). By looking at the wave frontthe speed is approximately 3.7m/s. This speed is radicalbig3µ/ρ which was alsoreported 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 ofthe phantom have shifted slightly over time.The same phantom was used for a second experiment. This time, insteadoftheplate, asmallscrewwasconnectedtothevoicecoil, whichactedalmostlike a point source. The result for this case is presented in Fig. 2.6(b). Thespeed from this figure is approximately 2.3m/s. This is the wave speedCatheline et al. reported in [36], i.e. the shear wave speed radicalbigµ/ρ.Note that in both cases Fig. 2.6(a) and (b), just after the compression isapplied, a compressional wave front is generated which propagates throughthe phantom at the speed of ultrasound, i.e. 1420m/s. The initial propaga-tion of this wave front cannot be tracked by the ultrasound, since the speedof ultrasound is the same as the wave speed. This wave front probes theboundaries of the medium and gets reflected back. However the displace-ments caused by multiple reflections of this wave from the boundaries are35Chapter 2. Elastography Using Longitudinal WavesFigure 2.6: (a) The result of the excitation with a plate which covers the topsurface 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 materialfor elastography studies, as its Young’s (and hence shear) modulus can beadjusted over a wide range (15–150kPa), without changing its density andthe speed of ultrasound cinf (Table 2.1). Figure 2.7 summarizes the Young’smodulus obtained vs. the volume portion of hardener or softener used inthe fabrication process for this material. However, the relaxation time andits frequency dependent behavior is not dependent on the volume portionof Softener or Hardener, as the graphs in Fig. 2.5(b) almost overlap for thefive test blocks. Therefore without any additives, this material by itself isunsuitable 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 mea-surements can be used to obtain all the elastic parameters. Other elasticparameters, namely Lam´e constants λ and µ, and Poisson’s ratio ν could beeasily found using the relations (2.7) to (2.9). The parameters are summa-rized in Table 2.6 for the five phantoms. The values of the Young’s modulusfrom the model fitting method were used in the calculations.The speed of ultrasound cinf has little variability among different softtissue types which are usually scanned using ultrasound imaging systems.36Chapter 2. Elastography Using Longitudinal WavesFigure 2.7: Obtained Young’s modulus as a function of the volume portion ofPlastic Softener or Hardener added during the fabrication of the phantomsTable 2.6: Summary of the elastic properties of the five test phantomsphantom λ µ E 0.5 - ν KNo (GPa) (kPa) (kPa) 1E-6× (GPa)1 2.1 6.0 18 1.5 2.12 2.1 8.7 26 2.1 2.13 2.0 18 55 4.7 2.04 2.1 26 77 6.1 2.15 2.0 44 132 11 2.037Chapter 2. Elastography Using Longitudinal WavesThis fact allows the conversion of time scale of ultrasound RF-pulses todistance for ultrasound image formation. This was also true for the SuperSoft PlasticTM used in our experiments. This fact has another importantconsequence which has been noted by other scholars as well [3]. Since softtissue is nearly incompressible, λ greatermuch µ; hence cinf is determined mainly by λ.In other words µ has no contribution. Since the variability of density is alsovery low (in living tissue and phantoms), λ should be almost constant. As aresult, the first Lam´e constant, λ, is not an elastic parameter of interest forcreating contrast in elastograms. In contrast, tissue stiffness is determinedby 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 elastographymethods based on shear wave excitation and methods based on longitudinalwave excitation. In shear wave excitation, shear wave speed cS = radicalbigµ/ρand hence shear modulus, µ is sought. In longitudinal excitation, Young’smodulus E = 3µ is sought. Hence these methods, although different inmany aspects, essentially use the same underlying physical quantity to createelastogram contrast.Another important point is that at low frequencies of excitation thephase difference in oscillations of the different points is low because of thelarge wavelength. Hence it is more difficult if not impossible to estimate theviscoelastic properties by using a low frequency. We conclude that a lowerlimit for the frequency of excitation exists. Note that this statement onlyapplies to the estimation of the absolute values of the viscoelastic propertiesusing the wave equation; relative estimations are still possible though. As amatter of fact, for a static compression ∂u/∂t ≡ ∂2u/∂t2 ≡ 0. In this casethe wave equation (2.29) implies ∂2u/∂z2 ≡ 0. The wave equation itself willbecome,ρ×0 = E ×0, (2.48)which cannot be inverted. However the strain, which is given by epsilon1 = ∂u/∂zis not necessarily equal to zero. Assuming a uniform stress, the strain isinversely proportional to the Young’s modulus E, and hence relative valuesof E can be obtained from the relative values of epsilon1. If in addition, forcemeasurement is carried out, the relative values can be turned into absolutevalues.It was observed that one-dimensional models would not fit the data verywell at higher frequencies of excitation. It is known that a block of materialacts as a wave-guide for elastic wave propagation [32], and the wave speed38Chapter 2. Elastography Using Longitudinal Wavesdepends on the dimensions of the source (the excited modes of the waveguide) as well as the frequency of excitation. Further study of these effectsis needed.Finally, the assumption of a 1D model which was made for the dis-placements is a limiting assumption. However, the methods presented inthis article are intended for rheometry purposes, i.e. characterization ofthe viscoelastic behavior of a test block assuming homogeneity. The tradi-tional force-displacement rheometry methods treat the entire test block asone unit, while in the presented methods an insight into internal vibrationsinside the block is obtained as well. Moreover no force measurements arecarried out. The presented methods do not yield local viscoelastic propertiesand hence are not presented as an alternative to any of the techniques in-tended for that purpose. Many research groups are looking at this problem,i.e. solving the inverse problem, from different perspectives. In particularthey use a more complete model, such as a 2D or 3D model, for the localinversion. This article, however, shows that longitudinal waves can also be avalid candidate for forming elastograms and extending the inverse problemmethods to this type of excitation.2.7 Summary and ConclusionIn this article an experimental framework was presented for viscoelasticcharacterization of tissue and tissue mimicking material. Longitudinal waveexcitation was used on five blocks of test phantoms and displacements wereimaged by an ultrasound system. Two methods were presented for analyz-ing the displacements. The first method estimated the Young’s modulusbased on the distance from the peak to the node of the standing wave pat-terns inside the block. The second method used a model fitting scheme tofind the values of the Young’s modulus and the relaxation time. The re-sults were compared to conventional rheometry results obtained using force-displacement measurements. For the elasticities good agreement was ob-served which validated the theory. For the viscosities, only the trends couldbe compared as the measurements were over different frequency ranges.We conclude that the speed of the longitudinal waves in finite media canbe orders of magnitude lower than their speed in infinite media and in factcomparable to the speed of shear waves. It is possible to track these wavesby using the ultrasound-based motion tracking techniques. The longitudinalwaves cannot simply be dismissed because of their presumed high speed ofpropagation or their insignificant amplitude. The viscoelastic properties39Chapter 2. Elastography Using Longitudinal Wavesof homogeneous phantom blocks can be estimated by applying the methodspresented in this article on the longitudinal displacements. The methods arebased on fitting a longitudinal-wave-equation-based model to the measureddisplacements. The estimated properties would be in agreement with thegold standard values from the conventional rheometry tests. Whether thehuman body or a particular tissue sample in vivo or ex vivo can be regardedas a finite medium, however, depends on the geometry of the setup and theexcitation used and needs further analysis.40References[1] A. Sarvazyan, “Handbook of elastic properties of solids, liquids andgases,” vol. III, 2001.[2] K. Parker, L. Gao, R. Lerner, and S. Levinson, “Techniques for elasticimaging: 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 imag-ing 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 unifiedview 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 ultra-sound,” Journal of Medical Ultrasonics, vol. 29, pp. 155–171, 2002.[6] R. Lerner and K. Parker, “Sonoelasticity images, ultrasonic tissue char-acterization and echographic imaging,” in Proceedings of the 7th Euro-pean Communities Workshop, 1987, nijmegen, Netherlands.[7] R. Lerner, K. Parker, J. Holen, R. Gramiak, and R. Waag, “Sono-elasticity: Medical elasticity images derived from ultrasound signals inmechanically vibrated targets,” Acoustical Imaging, vol. 16, pp. 317–327, 1988.[8] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastog-raphy: A quantitative method for imaging of elasticity of biologicaltissues,” Ultrason. Imaging, vol. 13, pp. 11–134, 1991.[9] O. Rudenko, A. Sarvazyan, and S. Emelianov, “Acoustic radiation forceand streaming induced by focused nonlinear ultrasound in a dissipativemedium,” J. Acoust. Soc. Am., vol. 99, no. 5, pp. 2791–2798, May 1996.41Chapter 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 medicaldiagnostics,” 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 modulusimaging with 2-d transient elastography,” IEEE Trans. Ultrason., Fer-roelect., Freq. Contr., vol. 49, no. 2, pp. 426–435, April 2002.[13] L. Sandrin, D. Cassereau, and M. Fink, “The role of the coupling termin transient elastography,” J. Acoust. Soc. Am., vol. 115, no. 1, pp.73–83, January 2004.[14] M. Fatemi and J. Greenleaf, “Ultrasound stimulated vibro-acousticspectrography,” Science, vol. 280, pp. 82–85, 1998.[15] K. Hoyt, K. J. Parker, and D. J. Rubens, “Real-time shear velocityimaging 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, andD. Bryk, “Acoustic properties of polyvinyl chloride gelatin for use inultrasonography,” Radiology, pp. 215–216, 1984.[17] T. J. Hall, M. Bilgen, M. F. Insana, and T. A. Krouskop, “Phantommaterials 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 ofpoly(vinyl alcohol) hydrogels produced by conventional crosslinking orby freezing/thawing methods,” Advances in Polymer Science, vol. 153,pp. 37–65, 2000.[19] E. Madsen, G. Frank, T. Krouskop, T. Varghese, F. Kallel, andJ. Ophir, “Tissue-mimicking oil-in-gelatin dispersions for use in het-erogeneous elastography phantoms,” Ultrasonic Imaging, vol. 25, pp.17–38, 2003.42Chapter 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, “Anthropomorphicbreast 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 dis-persions 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, “Acousticalproperties 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. Mon-grain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel me-chanical properties with four ultrasound elastography methods andcomparison with gold standard testings,” IEEE Trans. Ultrason., Fer-roelect., Freq. Contr., vol. 54, no. 9, pp. 498–509, March 2007.[24] U. Hamhaber, F. Grieshaber, J. Nagel, and U. Klose, “Comparison ofquantitative shear wave mrelastography with mechanical compressiontests,” Magn. Reson. Med., vol. 49, pp. 71–77, 2003.[25] H. Mansy, J. Grahe, and R. Sandler, “Elastic properties of syntheticmaterials 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. Mon-grain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel me-chanical properties with four ultrasound elastography methods andcomparison with gold standard testing,” IEEE Trans. Ultrason., Fer-roelect., Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007.[27] M. Doyley, J. Weaver, and K. Paulsen, “Measuring the mechanicalproperties of soft tissue over a wide frequency range by employing theprinciple of time temperature super-position,” in Proceedings of theSeventh International Conference on the Ultrasonic Measurement andImaging of Tissue Elasticity, Oct 2008, p. 75.[28] S. Srinivasan, T. Krouskop, and J. Ophir, “A quantitative comparisonof modulus images obtained using nanoindentation with strain elas-tograms,” Ultrasound Med. Biol., vol. 30, no. 7, pp. 899–918, 2004.43Chapter 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 androtary shear deformations,” Stud. Health Technol. Inform., vol. 94, pp.37–43, 2003.[30] V. Egorov, S. Tsyuryupa, S. Kanilo, M. Kogit, and A. Sarvazyan, “Softtissue elastometer,” Med. Eng. Phys., vol. 30, pp. 206–212, 2008.[31] A. Samani, J. Bishop, C. Luginbuhl, and D. B. Plewes, “Measuringthe elastic modulus of ex vivo small tissue samples,” Phys. Med. Biol.,vol. 48, pp. 2183–2198, 2003.[32] H. Kolsky, Stress Waves in Solids. New York: Dover Publications,Inc., 1963.[33] K. F. Graff, Wave Motion in Elastic Solids. Ely House, London:Oxford University Press, 1975.[34] F. Lockett, Non-linear Viscoelastic Solids. London: Academic Press,1972.[35] W. Flugge, Viscoelasticity, 2nd ed. Berlin: Springer, 1975.[36] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelka-ram, and J. Culiolic, “Measurement of viscoelastic properties of ho-mogeneous soft solid using transient elastography: An inverse problemapproach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec2004.[37] F. A. Duck, Physical Properties of Tissue: A Comprehensive ReferenceBook. Academic Press, 1990.[38] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complex-valued stiffness reconstruction by for magnetic resonance elastographyby algebraic inversion of the differential equation,” Magn. Reson. Med.,vol. 45, pp. 299–310, 2001.[39] R. S. C. Cobbold, Foundations of Biomedical Ultrasound. OxfordUniversity Press, 2006.[40] H. Eskandari, S. E. Salcudean, and R. Rohling, “Viscoelastic parameterestimation based on spectral analysis,” IEEE Trans. Ultrason., Ferro-elect., Freq. Contr., vol. 55, no. 7, pp. 1611–1625, 2008.44Chapter 2. Elastography Using Longitudinal Waves[41] R. Zahiri and S. Salcudean, “Motion estimation in ultrasound imagesusing time-domain cross-correlation with prior estimates,” IEEE Trans.Biomed. Eng., vol. 53, no. 10, pp. 1990–2000, 2006.45Chapter 3Theoretical Limitations ofthe Elastic Wave EquationInversion for TissueElastography∗3.1 IntroductionElastography is a fast developing field of medical imaging which strives toprovide images of the mechanical properties of tissue [1–8]. The main me-chanical property of interest is the elasticity or Young’s modulus E. Othermechanical properties of interest include viscosities, relaxation times, non-linearity parameters, harmonics, poro-elasticity and Poisson’s ratio [8–22].The majority of the elastography techniques have three components incommon:• an excitation mechanism which creates quasi-static or dynamic defor-mations in the tissue [23–31],• a medical imaging system augmented by a motion estimation tech-nique which is capable of providing displacement (or velocity) imagesof the tissue while it is being deformed [25,32–41], and• an inversion technique which transforms the displacement images intoelasticity 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 typicallyuse 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. 200946Chapter 3. Theoretical Limitationsfind the best local elasticity values which would create the same dis-placement pattern as those observed inside the tissue, under similarboundary and excitation conditions. To find the optimum elasticityvalues, the forward problem is solved iteratively. In each iterationthe elasticity values are adjusted until the computed displacementsfrom the model match the measured displacement from the tissue. Inaddition to the measured displacements, these techniques typically re-quire further information about the shape of the tissue, its boundaryconditions, and the excitation.• Local inversion techniques [11,25,27–30,37,50–52,58–61]: Typicallyin these techniques, a particular wave equation is assumed to governthe displacements. This partial differential equation (PDE) relatesthe spatial derivatives to the temporal derivatives of the displacementwith the local elasticity as the coefficient. If the displacement mea-surements are available over a range of space and time, the elasticitiescan be found either by an algebraic inversion of the wave equation orby estimating the phase speed from the gradient of the phase. Thenext 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 forestimating the tissue elasticity.• Study the effects of the formulation of the wave equation on the in-version 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 con-tinuum is studied. The hypothesis to be proven is that the phase speed,which is the final estimated quantity in many of the local inversion tech-niques, does not exactly represent the local mechanical properties of thetissue. The hypothesis is proven by showing that in addition to the localmechanical properties of the tissue, the phase speed depends on the geom-etry of the wave as well as the geometry of the tissue itself. This effect isknown as the Lamb wave effect specially in the context of thin plates [62].The practical implication is that artifacts are inevitably present in the elas-tograms, no matter how well the algorithm is implemented or how accuratethe displacement measurements are. Nevertheless, many of the elastographytechniques have been proven clinically to provide valuable information aboutthe elasticity of the tissue. The excitation in these techniques is chosen sothat the actual displacements created in the tissue satisfy the assumptionsmade, which in turn justifies the inversion techniques used. These meth-47Chapter 3. Theoretical Limitationsods include quasi-static constant strain compression, pulsed excitation withexternal exciters, and acoustic-radiation-force-based techniques.A previous study of the limitations of the local inversion techniques forestimating the tissue elasticity can be found in [63]. In that work the au-thors experimentally measured the displacements caused by a circular pistoninside 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 regimesof vibration, a region was identified when transient shear waves are prop-agating in the medium and phase speed measurements can be carried out.The deleterious effects of the medium boundaries were dealt with by reject-ing the steady state regime for phase speed measurements. The effects of theboundary conditions imposed by the exciter were also studied. An approxi-mate model, the Rayleigh–Sommerfeld’s solution, and the Green’s functionwere used in this analysis. They concluded that the size of the exciter af-fects the measured phase speed. Another conclusion was that at very lowfrequencies, the effects of the longitudinal wave on the phase speed cannotbe neglected. The final preferred method was to use a pulsed transient exci-tation with a small piston to get more accurate viscoelastic measurements.The analysis of the aforementioned article is limited to the particularconfiguration used, which is very similar to a point source on a semi-infinitemedium. In many elastography applications, to get a desired level of dis-placements inside tissue, it is of interest to use exciters with a larger foot-print. We have shown previously that if the size of the exciter is comparableto the size of the sample, the phase speed observed is the extensional wavespeed [64]. In the present article we look at the problems associated withphase speed measurements from a more general point of view, and discusssome of the phenomena that affect the accuracy of the measurements. Wewill show that neglecting these phenomena results in un-accounted-for arti-facts in the computed elastograms.We start by a brief review of the state-of-the-art local inversion tech-niques in Section 3.2. Section 3.3 is a critique on the use of the terminologyassociated with the elastic waves. The main theoretical results of the articleare presented in Section 3.4. Simulation results are presented in Section 3.5.Section 3.6 is devoted to the conclusion.48Chapter 3. Theoretical Limitations3.2 Local Inversion of the Wave EquationThe mechanical wave equation in an isotropic elastic continuum written inthe Cartesian coordinate system is [62]ρ∂2ui∂t2 = (λ+µ)∂∂i∆+µ∇2ui i = x,y,z , (3.1)where x = [x,y,z]T are the Cartesian coordinates which define the spatiallocation of the point, t is the time variable, u(x,t) = [ux,uy,uz]T is thedisplacement field, and the coefficients ρ(x), λ(x) and µ(x) are the density,the first, and the second Lam´e constants respectively. Here, we have adoptedthe notation of Kolsky [62]. The dilatation ∆ is defined by∆ = ∂ux∂x + ∂uy∂y + ∂uz∂z , (3.2)and determines the volumetric changes as the wave propagates. ∇2 is theLaplacian operator:∇2 = ∂2∂x2 +∂2∂y2 +∂2∂z2 . (3.3)We denote the Poisson’s ratio and the Young’s modulus by ν and E respec-tively.The material of interest in the field of elastography is tissue in vivo. Thismaterial is found to be nearly incompressible [51,65], hence the values of∆ are typically very small. Other implications of the incompressibility interms of the mechanical properties are as follows:ν ≈ 0.5 , λ greatermuch µ 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 shouldbe based on a physical quantity which has a high degree of variability amongdifferent tissue types. The density of most soft tissue is close to the densityof water [66], 1000kg/m3 so ρ is not a candidate. Neither is the first Lam´econstant λ which is found to be close to the λ of water [67], i.e. 2.3GPa. Onthe 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 elasticityE = 3µ. Given the displacements as a function of time and space, the49Chapter 3. Theoretical Limitationsmost that can be recovered from equation (3.1) are the ratios λ/ρ and µ/ρ.In practice, the elastograms are based on E/ρ. Since the variability of ρ issmall, the obtained contrast is almost the same as the one obtained by usingthe absolute value of elasticity E.A number of issues arise when medical imaging devices are used to es-timate the displacement fields. In many cases it is not possible to acquireall three components of the displacement field, i.e. ux, uy, and uz. In othercases the acquired displacement field may not be available over a volume.For instance, in most cases of ultrasound elastography with 2D probes, onlyux(x,y) is available and only on a single plane z =constant and not overa volume. Some customized ultrasound systems can estimate the displace-ments over a volume and in multiple directions [37,41,72]. MRI basedtechniques on the other hand, can acquire the full displacement field over avolume, at the expense of longer acquisition times [14,25,50,51,73–75].Different approaches to the inversion of the wave equation have beentaken. As mentioned earlier the dilatation ∆(x,t) levels in soft tissue arewell beyond the accuracy of the measurement systems. Therefore the waveequation (3.1) cannot directly be inverted to obtain µ/ρ. Some state-of-the-art 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ρ∂2ui∂t2 = µ∇2ui i = x,y,z. (3.5)Note that for an anisotropic inversion all three components are nec-essary. If isotropy is assumed, measuring one component of the dis-placement field, for instance ux, is enough. In this case, it is naturalto chose the measurement axis (the x-axis) aligned with the axis onwhich the wave causes the greatest displacements. If the measurementis only available on a plane, for instance ux(x,y), or on a single line,50Chapter 3. Theoretical Limitationsux(x), then the following assumptions are usually made∇2ux =parenleftbigg ∂2∂x2 +∂2∂y2 +∂2∂z2parenrightbiggux≈parenleftbigg ∂2∂x2 +∂2∂y2parenrightbiggux≈ ∂2∂x2ux . (3.6)• Assuming a shear wave equation [52]: In this case it is assumed thatthe wave is a plane shear wave. If the wave causes the greatest dis-placements in the x-direction, then this component of the displace-ment field is measured and assumed to satisfy a 1D shear wave equa-tion,ρ∂2ux∂t2 = µ∂2ux∂y2 . (3.7)• Assuming a thin rod (extensional) wave equation [64]: When the sizeof a compressing exciter is comparable to the dimensions of the tis-sue being imaged the tissue behaves like a thin rod experiencing ex-tensional waves. Under these conditions the component of the dis-placement field parallel to the direction of the excitation, say the x-component, satisfiesρ∂2ux∂t2 = E∂2ux∂x2 . (3.8)• Writing the shear wave equation for the rotations [14,53,72–75]: De-fine the rotation w(x,t) as half the curl of the displacement field:w(x,t) = 12∇×u(x,t) (3.9)From (3.1) each component of the rotation field satisfies a separateshear wave equation,ρ∂2wi∂t2 = µ∂2wi∂i2 i = x,y,z. (3.10)Note that this approach does not make any assumptions, but it re-quires the measurement over a 3D volume of at least two componentsof the displacement field for an isotropic inversion. For an anisotropicinversion, all the components are necessary.The local inversion algorithms (with the exception of the last approach)51Chapter 3. Theoretical Limitationsreduce the general form of the wave equation (3.1) to a 1D wave equationof the formρ∂2ux∂t2 = f(E)∂2ux∂x2 , (3.11)where f(E) is a function of the local elasticity. The quantitycph =radicalbigf(E)/ρ (3.12)in this 1D equation is called the phase speed. For a unidirectional harmonicsolution,ux(x,t) = aexp(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)ρ =∂2ux/∂t2∂2ux/∂x2 , (3.14)This method is sometimes called the Algebraic Inversion of the DifferentialEquation (AIDE) [50]. Another method is to use the concept of the phasespeed, radicalbigf(E)/ρ = cph = ω∂∂x∠ux, (3.15)sometimes called the Phase Gradient method [50].If f(E) were only a function (even an unknown function) of the localmechanical properties of the tissue, the elastograms based on the phasespeed,radicalbigf(E)/ρ would be free from artifacts. However as we show in Section3.4, f(E) is not only a function of the local mechanical properties of themedium but also a function of the geometry of the wave and the geometryof 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 phasespeed can coexist.3.3 A Critique on the Terminology Associatedwith the Wave EquationBefore we move on to our results on the phase speed we present our choiceof the terminology. By the fundamental theorem of vector calculus, orHelmholtz decomposition theorem, any sufficiently smooth and decaying52Chapter 3. Theoretical Limitationsvector field can be written as the sum of a divergence-free and a curl-freevector field [76]. The solutions of the wave equation (3.1) satisfy the condi-tions of the theorem, and therefore it is possible to write a vibration patterninside an elastic medium as the sum of a divergence-free and a curl-freedisplacement field.• Dilatational wave: or better termed irrotational component is a com-ponent of the wave which is curl-free. This component of the wavephenomena is solely the result of voluminal changes in the medium.• Distortion wave: or better termed equivoluminal component is a com-ponent of the wave which is divergence-free. This component of thewave phenomena preserves the volumes of the infinitesimal elements.There is some confusion however, over the definition and applicationof 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 theparticle velocities are parallel to the phase velocity. For a harmonicwave the particle velocity would be perpendicular to the equiphasesurfaces.• Transverse wave: or shear wave is a wave for which the particle ve-locities are perpendicular to the phase velocity. For a harmonic wavethe particle velocity would be tangent to the equiphase surfaces.These definitions are very intuitive with respect to the literal meanings of theterms. However a major problem exists with these definitions; to the bestknowledge of the authors, there is no theorem proving that the quality of awave being transversal or longitudinal (taken these definitions) is invariantwith respect to the wave equation. In other words, a wave propagating in aninfinite elastic medium might be purely longitudinal at one instant but notso at a later instant (even without encountering a boundary and undergoingmode conversion). This severely limits the applicability of these terms.When these definitions are taken, a connection is usually assumed be-tween longitudinal and irrotational waves on one hand and transverse andequivoluminal waves on the other hand. The origin of this connection couldbe the study of the simple case of a plane wave in an infinite medium. In-deed a longitudinal plane wave propagating in an infinite medium has noequivoluminal components and consists solely of a dilatational component.This wave propagates at a phase speed of radicalbig(λ+2µ)/ρ hence the name lon-gitudinal wave speed. Similarly a shear plane wave propagating in an infinitemedium has no dilatational components and consists solely of an equivolu-minal component. This wave propagates at a phase speed of radicalbigµ/ρ hence53Chapter 3. Theoretical Limitationsthe name shear wave speed.A better way to define these terms, however, is to take for them the samedefinitions as presented for the “dilatational waves” and “distortion waves”.This is the accepted nomenclature in the physics literature, specially in thefield of the electromagnetic waves [77]:• Longitudinal wave: or compressional wave is a component of the wavewhich is curl-free.• Transverse wave: or shear wave is a component of the wave which isdivergence-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 thata wave can contain both the longitudinal and transversal components, ormay lack one. In some cases one is not interested in the decomposition ofthe displacement field into these components, but rather in the study of thetotal displacement itself. However since each of these components separatelysatisfies a reduced wave equation with a definite wave speed, it is at othertimes useful to study them separately. The speed of the longitudinal com-ponent and the transverse component is equal to radicalbig(λ+2µ)/ρ and radicalbigµ/ρrespectively.Since the wave equation can be reduced to the two aforementioned waveequations, it is evident that in an infinite elastic medium, a purely longi-tudinal wave will remain purely longitudinal, and a purely transverse wavewill remain purely transverse. This invariance property makes the termswell-defined. The drawback is that the intuitive relationship between thedirections 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 transversalcould be determined solely from the displacement vector field at a singleinstant. The wave equation preserves this quality as the time evolves.As will become clear in the next section, the relationship between thephase speed and the longitudinal and shear wave speeds, radicalbig(λ+2µ)/ρ andradicalbigµ/ρ, is a complex relationship which depends on the geometry of the wave,as well as the geometry of the medium.54Chapter 3. Theoretical Limitations3.4 Phase Speed of Mechanical Waves in ElasticMediums3.4.1 Dependence of the Phase Speed on the Geometry ofthe Wave: Infinite MediumsTheorem. Consider an infinite homogeneous linear elastic medium withdensity ρ and Lam´e constants λ and µ. Given the angular frequency ω fora harmonic wave:• for any number c such thatradicalbiggµρ ≤ c <radicalBiggλ+2µρ , (3.16)there exists a shear beam form for the harmonic wave for which thephase speed is equal to c.• for any number c such thatradicalBiggλ+2µρ ≤ c, (3.17)there exists infinitely many beam forms for the harmonic wave forwhich the phase speed is equal to c. These waves contain both longi-tudinal and shear components.The theoretical derivation which follows can be found in classical text-books on elastic waves [62,78] as part of the solution to the wave equation incylindrical coordinate systems. However studying the implications in termsof the phase speed, as summarized in the theorem, is novel.Proof. Consider the elastic wave equation in the cylindrical coordinate sys-tem,ρ∂2ur∂t2 = (λ+2µ)∂∆∂r −2µr∂¯ωz∂θ +2µ∂¯ωθ∂z , (3.18)ρ∂2uθ∂t2 = (λ+2µ)1r∂∆∂θ −2µ∂¯ωr∂z +2µ∂¯ωz∂r , (3.19)ρ∂2uz∂t2 = (λ+2µ)∂∆∂z −2µr∂∂r (r¯ωθ)+2µr∂¯ωr∂θ , (3.20)where (r,θ,z) are the cylindrical coordinates and u = (ur,uθ,uz) is the55Chapter 3. Theoretical Limitationsdisplacement field. The dilatation in the cylindrical coordinates is given by,∆ = 1r ∂∂r (rur)+ 1r ∂uθ∂θ + ∂uz∂z , (3.21)and the rotations are given by,2¯ωr = 1r ∂uz∂θ − ∂uθ∂z , (3.22)2¯ωθ = ∂ur∂z − ∂uz∂r , (3.23)2¯ωz = 1rbracketleftbigg ∂∂r (ruθ)−∂ur∂θbracketrightbigg. (3.24)Consider a cylindrical wave beam propagating along the z-axis with thecenter 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 assumethat the wave is symmetrical around the z-axis, hence ∂/∂θ annihilates thevariables. From (3.22) and (3.24) ¯ωr = ¯ωz = 0. The reduced wave equationin this case would be,ρ∂2ur∂t2 = (λ+2µ)∂∆∂r +2µ∂¯ωθ∂z , (3.25)ρ∂2uz∂t2 = (λ+2µ)∂∆∂z −2µr∂∂r (r¯ωθ) . (3.26)We are interested in harmonic waves propagating along the z-axis, forinstance in the negative z-direction. The general form of such a wave is,ur = U(r)exp(i(kzz +ωt)) , (3.27)uz = W(r)exp(i(kzz +ω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 isequal to cph = ω/kz. Substitution in (3.21) and (3.23) yields the expressions56Chapter 3. Theoretical Limitationsfor the dilatation and rotation:∆(r,z,t) =bracketleftbigg∂U(r)∂r +U(r)r +ikzW(r)bracketrightbiggexp(i(kzz +ωt)) , (3.29)2¯ωθ(r,z,t) =bracketleftbiggikzU(r)− ∂W(r)∂rbracketrightbiggexp(i(kzz +ωt)) . (3.30)The forms (3.27) and (3.28) simplify the wave equation (3.25) and (3.26),−ρω2ur = (λ+2µ)∂∆∂r +2iµkz¯ωθ , (3.31)−ρω2uz = ikz(λ+2µ)∆− 2µr ∂∂r (r¯ωθ) . (3.32)By eliminating ¯ωθ and ∆ between these equations, the longitudinal andtransversal wave equations are derived:∂2∆∂r2 +1r∂∆∂r +k∆2∆ = 0, (3.33)∂2¯ωθ∂r2 +1r∂¯ωθ∂r −¯ωθr2 +k¯ωθ2¯ωθ = 0, (3.34)where k∆ and k¯ωθ are the geometrical beam numbers defined byk∆2 = ω2λ+2µρ−kz2 , (3.35)k¯ωθ2 = ω2µρ−kz2 . (3.36)In the case that the desired phase speed, c, satisfies (3.16) the corre-sponding 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 variablesfrom r to k∆r and k¯ωθr respectively. The resulting equations are Besselequations of the zeroth and first orders respectively. The physically mean-57Chapter 3. Theoretical Limitationsingful 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 andfirst orders respectively. Now the forms for the dilatation and rotation aregiven 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,U(r) = C1 ∂∂rJ0(k∆r)+C2kzJ1(k¯ωθr)= −C1k∆J1(k∆r)+C2kzJ1(k¯ωθr), (3.40)W(r) = C1ikzJ0(k∆r)+C2 ir ∂∂r [rJ1(k¯ωθr)]= C1ikzJ0(k∆r)+C2ik¯ωθJ0(k¯ωθr), (3.41)where C1 and C2 are two arbitrary constants. From (3.29) and (3.30) thedilatation and rotation become∆(r,z,t) =−2C1(kz2 +k∆2)J0(k∆r)exp(i(kzz +ωt)) , (3.42)2¯ωθ(r,z,t) =2iC2(kz2 +k¯ωθ2)J1(k¯ωθr)exp(i(kzz +ωt)) . (3.43)The wave can thus be a mixture of the longitudinal and shear componentswith C1 and C2 defining the respective proportions.Now if k∆ is imaginary, the Bessel functions J0(k∆r) and J1(k∆r) go toinfinity as r goes to infinity. This is not physically meaningful. Thus for aphase 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 travelsat the phase speed c.On the other hand if the phase speed c satisfies (3.17), both C1 and C2can have nonzero values. Therefore infinitely many beam forms exist forwhich the wave travels at the phase speed c. Each of these beam formscontains both the longitudinal and shear components. This completes theproof.It follows that the phase speed depends not only on the mechanical58Chapter 3. Theoretical Limitationsproperties of the medium, but also on the geometry of the wave. Even fora purely shear wave (C1 = 0), the phase speed can have any value whichis greater than or equal to radicalbigµ/ρ. The shear wave speed, radicalbigµ/ρ, is not thephase speed of every shear wave in an infinite medium. It is however thephase speed of the uniform plane shear waves. This issue is in addition to themain drawback of the use of the phase speed [63] for tissue characterization;namely that the phase speed can not be defined when multiple waves aretraveling in different directions, for instance when there are reflections ofthe wave from the boundaries.3.4.2 Dependence of the Phase Speed on the Geometry ofthe Medium: Wave GuidesIn the previous section the medium was assumed to be infinite in size. Thisassumption 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 wepresent some classical results on the wave guides and study their implicationsin the field of elastography. The wave guides are infinite in at least onedirection and finite in at least another. Therefore they cannot model thetissue behavior accurately either. However they can be considered as anintermediate step in moving from the analysis of an infinite medium to abounded medium. As such the insight gained from studying them is usefulin understanding and designing elastography systems.The wave guide we study is an infinitely long cylindrical rod. Choosethe cylindrical coordinate system with the axis of the cylinder on the z-axis.As in the previous section, we are interested in harmonic waves propagatingalong the z-axis, i.e. along the axis of the cylinder, which are symmetricalaround 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) satisfiesthe wave equation with phase speed c, provided that kz is chosen to satisfyω/kz = c. However in this case the boundary conditions impose restrictionson the permissible values of c.The boundary of the cylinder is free from stresses. The expressions for59Chapter 3. Theoretical Limitationsstresses in the cylindrical coordinate system are,σrr = λ∆+2µ∂ur∂r , (3.44)σrθ = µbracketleftbigg1r∂ur∂θ +r∂∂rparenleftBiguθrparenrightBigbracketrightbigg, (3.45)σrz = µbracketleftbigg∂ur∂z +∂uz∂rbracketrightbigg. (3.46)By the assumptions σrθ vanishes everywhere. The boundary conditionsrequire that σrr and σrz vanish on the surface of the cylinder. Denotethe radius of the cylinder by a. Substitution of the expressions for ur anduz from (3.27), (3.28), (3.40), and (3.41) translates the boundary conditionsinto the following two equations,C1bracketleftbigg2µ ∂2∂r2vextendsinglevextendsinglevextendsinglevextendsingler=aJ0(k∆r)− λρω2λ+2µJ0(k∆a)bracketrightbigg+2C2µkz ∂∂rvextendsinglevextendsinglevextendsinglevextendsingler=aJ1(k¯ωθr) = 0, (3.47)2C1kz ∂∂rvextendsinglevextendsinglevextendsinglevextendsingler=aJ0(k∆r)+C2parenleftbigg2kz2 − ρω2µparenrightbiggJ1(k¯ωθa) = 0. (3.48)Note that these equations depend only on the ratio of C1/C2. Eliminatingthis variable between the two, a single equation is obtained for the unknownkz. This equation is called the Pochhammer frequency equation [62]. For agiven material ρ, µ, λ and a are known. So the equation basically describesthe relationship between the wave number kz and angular frequency ω ofthe wave, or equivalently the relationship between the phase speed c andω. As it turns out, unlike the infinite medium, not every phase speed ispossible for a given ω. The permissible speeds must satisfy the Pochhammerfrequency equation. Substitution into either of the equations (3.47) or (3.48)determines the ratio C1/C2, i.e. the proportions of the longitudinal andtransversal waves that should be mixed together to obtain that particularphase speed.60Chapter 3. Theoretical LimitationsTable 3.1: Mechanical property values used for the simulationρ λ µkg/m3 GPa kPa1000 2.3 103.5 Simulation Results3.5.1 Shear Waves in a Tissue Mimicking Infinite MediumTo gain some insight into the theoretical results which were presented inthe previous section, we simulate the shear wave beam forms in an infinitemedium for different angular frequencies and phase speeds. For a chosenpair of angular frequency ω and phase speed c, the wave number kz is foundfrom (3.37). We choose the mechanical properties of the medium to matchthose found in the soft living tissue such as the breast. These values arelisted in Table 3.1. Given the mechanical properties of the medium, k∆ andk¯ωθ are found from (3.35) and (3.36) respectively. In the next step U(r) andW(r) are found from (3.40) and (3.41), by setting C1 = 0 and choosing anarbitrary value for C2. Note that C1 is chosen to be zero to obtain a purelyshear wave and the value of C2 only determines the overall amplitude andphase of the wave, but does not change the waveform. If the time evolutionof 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).Asshearwavescanpropagateatanyphasespeedaboveradicalbigµ/ρ = √10m/s,the choice of the phase speed, c, becomes arbitrary. We simulated the beamforms using two chosen values of c = 5m/s and c = 12m/s. The resultsare shown in Fig. 3.1 and 3.2 respectively. The first and the second col-umn in these figures show the displacement components ur((r,θ,z),t) anduz((r,θ,z),t) respectively at t = 0. The three rows correspond to increasingfrequencies 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 sameharmonic motion as ur((r,θ,z),t) and uz((r,θ,z),t) over an infinite crosssection of the medium, for instance at z = 0, the wave forms due to thiscylindrically symmetric infinite exciter, in the steady state, would be thesame as those depicted in these figures. The figures show how changingthe excitation pattern of such an exciter can result in a completely differentphase speed, even if the frequency of the excitation is not changed.61Chapter 3. Theoretical LimitationsFigure 3.1: Infinite medium; shear beam patterns for three different frequen-cies all sharing the same phase speed of 5m/s. (a), (b), (c) radial componentsof the displacement (d), (e), (f) z-component of the displacement at 40Hz,65Hz, and 100Hz respectively. Note that only a cubic portion of the mediumnear the axis of symmetry is shown; The medium extends in all directionsto infinity.62Chapter 3. Theoretical LimitationsFigure 3.2: Infinite medium; shear beam patterns for three different fre-quencies all sharing the same phase speed of 12m/s. (a), (b), (c) radialcomponents of the displacement (d), (e), (f) z-component of the displace-ment at 40Hz, 65Hz, and 100Hz respectively. Note that only a cubic portionof the medium near the axis of symmetry is shown; The medium extends inall directions to infinity.63Chapter 3. Theoretical LimitationsTable 3.2: Mechanical property values used for the simulationρ λ µ akg/m3 GPa kPa m1000 2.3 10 0.053.5.2 Waves in a Tissue Mimicking Cylindrical RodThe Pochhammer equation has been studied for metallic rods such as steelbeams. To gain some insight into this equation when dealing with the liv-ing tissue, we consider a simple example here. We choose the mechanicalproperties of the medium to match those found in the living tissue. Thesevalues are listed in Table 3.2. Using these values, the Pochhammer equationwas solved numerically using MATLAB for different frequencies. Figure 3.3shows the plots of the values of c obtained for each ω. At higher frequen-cies, the equation has multiple roots (modes), thus the multiples plots inthis figure.Mode 0 has a constant phase speed c = radicalbigµ/ρ, i.e. the shear wave speedfor 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 trivialsolution. As a matter of fact assuming a phase speed equal to the shearwave speed results in vanishing displacements, independent of the materialproperties. The implication is that no (axis symmetric) wave can travelalong the cylindrical rod with the shear wave speed.Figure 3.3: The first four modes of the cylindrical bar64Chapter 3. Theoretical LimitationsAt 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 isin contrast to the case of the infinite medium studied before, in which thewaves, independent of their frequency, could travel at a range of speeds. Thelow frequency speed is equal to radicalbigE/ρ which is the familiar value obtainedfrom the thin rod approximation theory [64]. This also shows the theoreticaljustification behind the assumptions made in (3.8). As the frequency goeshigher, other modes start to appear.At high frequencies (above 38Hz) multiple waves of the same frequencycan propagate simultaneously, each with a different phase speed. For in-stance 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 inFig. 3.4. In this case it may not be possible to recover a phase speed fromstudying the displacement patterns inside the material.To see this more clearly, assume that a harmonic exciter vibrating ata frequency of 100Hz is placed at infinity on the rod and has caused thefollowing wave pattern:ur(t) = 3Uc1(r)exp(i(207.1z +2pi100t))+5Uc2(r)exp(i(168.94z +2pi100t))+8Uc3(r)exp(i(111.60z +2pi100t)) , (3.49)uz(t) = 3Wc1(r)exp(i(207.1z +2pi100t))+5Wc2(r)exp(i(168.94z +2pi100t))+8Wc3(r)exp(i(111.60z +2pi100t)) , (3.50)where the Uci and Wci are the solutions presented in Fig. 3.4 for the threephase speeds. Because of the linearity of the wave equation, any linearcombination of the solutions presented in Fig. 3.4 is a solution. In partic-ular the above linear combination with the coefficients 3, 5, and 8 causedby the particular exciter used, is a solution. It is not hard to verify bysubstitution that neither the direct inversion method (3.14), nor the phasegradient method (3.15) results in a meaningful mechanical property for thehomogeneous cylinder.The implication in the context of elastography is that at higher frequen-cies, multiple modes appear in the measured displacements which makes itimpossible to recover a single phase speed and determine the mechanicalproperties based on that.65Chapter 3. Theoretical LimitationsFigure 3.4: Permissible beam patterns (modes) in an infinite cylindrical rodfor waves having a frequency of 100Hz. Three (and only three) phase speedsare 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.66Chapter 3. Theoretical Limitations3.6 ConclusionThe propagation of elastic waves trapped in finite media, such as soft tissue,is a complex phenomenon even without considering the viscous and non-linear effects. The term wave is generally applied to any vibration patterninside the tissue. However the phase velocity cannot be defined for all thecases of the vibration patterns. Usually multiple waves traveling in differ-ent directions superimpose to create the vibration pattern. In this case thephase velocity might be well-defined for each wave, but not for the resultantof the superposition.Even in the simple case where a single frequency wave is traveling in asingle direction, the geometry of the wave and the medium create multiplepermissible phase speeds, and thus make it difficult to recover a meaningfulphase speed for the traveling wave. The parameters used in our simulationswere chosen to be close to those of the living tissue. The diameter of thecylindrical rod was chosen to be 10cm which is in the same size range ashuman organs. Yet multiple phase speeds and modal behavior start toappear at frequencies as low as 38Hz. This frequency is in the range offrequencies used in many dynamic excitation elastography techniques.One remedy seems to be the use of the natural decomposition of thevibrations into the longitudinal (dilatational) and transverse (shear) com-ponents. Each of these components always satisfies its own wave equationwith its own wave speed which is independent of the geometry (note thatthis is not the phase speed however). Taking the divergence of the measureddisplacement field yields the dilatation which satisfies the dilatational waveequation with wave speed radicalbig(λ+2µ)/ρ,ρ∂2∆∂t2 = (λ+2µ)∇2∆. (3.51)Since living tissue is incompressible, the dilatations will be very minuteand impossible to detect by taking the spatial derivatives of the measureddisplacements from any of the currently available imaging techniques. Evenif accurate measurements were possible, inverting this wave equation wouldresult in the local values of the longitudinal wave speed, radicalbig(λ+2µ)/ρ, beingknown. 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 theformed elastogram would be minimal.On the other hand, taking the curl of the displacement field results in67Chapter 3. Theoretical Limitationsthe rotation fields,∇×(ux,uy,uz)T = (¯ωx, ¯ωy, ¯ωz)T (3.52)each of which satisfies a shear wave equation with wave speed radicalbigµ/ρ,ρ∂2¯ωi∂t2 = µ∇2¯ωi i = x,y,z. (3.53)Since living tissue is incompressible, E ≈ 3µ and inverting these equationsresults in local information for tissue elasticity E. Moreover since theseequations are naturally 1D, it is possible to use the phase speed in theircontext.From this discussion it is concluded that the only theoretically flawlessmethod of imaging tissue elasticity is by taking the curl of the displacementfields. This requires the measurement of the displacements in all three direc-tions over the volume of interest. In cases where this is not feasible, such asultrasound elastography, artifacts will always be present. These artifacts arenot just due to noisy measurements, poor algorithms, or imperfect setupsbut 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 directionat a single speed:• The excitation amplitude should be chosen small enough to reduce thereflected wave amplitudes to the level of other measurement noises.The damping present in the tissue helps by reducing the amplitude ofthe reflected waves from the boundaries of the tissue and bones.• The frequency of excitation should be kept to a minimum to preventthe appearance of higher modes. However lower frequencies result inlower damping, and therefore a trade off exists here.• The excitation pattern should be chosen so as to excite mainly thelowest mode.From an engineering point of view, however, many factors beyond theanalysis presented here affect the success of an elastography system. Many ofthe devised elastography techniques are proven to provide clinically valuableimages and some of them have even found their way into the market. Quasi-static constant strain, pulsed excitation with external exciters, acoustic-radiation-force-based and other techniques are each designed to make thedisplacements created in the clinical setting match the assumptions madeabout them. This results in the validity of their inversion techniques and68Chapter 3. Theoretical Limitationstheir valuable elastograms.The presence of the artifacts, by itself, does not mean that an elastog-raphy system should be dismissed. After all, each of the available medicalimaging modalities has its own artifacts, and yet they are very well provento be useful in diagnosis and treatment.69References[1] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastog-raphy: A quantitative method for imaging of elasticity of biologicaltissues,” Ultrasonic Imaging, vol. 13, pp. 111–134, 1991.[2] K. Parker, L. Gao, R. Lerner, and S. Levinson, “Techniques for elas-tic imaging: A review,” IEEE Engineering in Medicine and BiologyMagazine, vol. 15, no. 6, pp. 52–59, 1996.[3] L. Gao, K. J. Parker, R. M. Lerner, and S. Levinson, “Imaging of theelastic 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 andGases. Academic Press, 2001, vol. III, ch. Elastic Properties of SoftTissues.[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 ultra-sound,” Journal of Medical Ultrasonics, vol. 29, pp. 155–171, 2002.[7] J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imag-ing elastic properties of biological tissues,” Annual Review BiomedicalEngineering, vol. 5, pp. 57–78, April 2003.[8] K. J. Parker, L. S. Taylor, S. Gracewski, and D. J. Rubens, “A unifiedview 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 re-lationships in tissue and their effect on the contrast-to-noise ratio inelastograms,” Ultrasound Med. Biol., vol. 26, no. 5, pp. 839–851, 2000.70Chapter 3. Theoretical Limitations[10] E. E. Konofagou, T. P. Harrigan, J. Ophir, and T. A. Krouskop, “Poroe-lastography: Imaging the poroelastic properties of tissues,” UltrasoundMed. Biol., vol. 27, no. 10, pp. 1387–1397, 2001.[11] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity andviscosity 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 inthe impulse diffraction field of elastic waves induced by the acoustic ra-diation 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 feasi-bility of using elastography for imaging the poisson’s ratio in porousmedia,” 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 mrelastography,” Magnetic Resonance Imaging, vol. 23, pp. 159–165, 2005.[15] R. Righetti, J. Ophir, and T. A. Krouskop, “A method for generat-ing permeability elastograms and poisson’s ratio time-constant elas-tograms,” Ultrasound Med. Biol., vol. 31, no. 6, pp. 803–816, 2005.[16] S. Chen, R. Kinnick, J. F. Greenleaf, and M. Fatemi, “Difference fre-quency and its harmonic emitted by microbubbles under dual frequencyexcitation,” Ultrasonics, vol. 44, pp. 123–126, 2006.[17] R. Sinkus, J. Bercoff, M. Tanter, J.-L. Gennisson, C. E. Khoury, V. Ser-vois, A. Tardivon, and M. Fink, “Nonlinear viscoelastic properties oftissue 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, “Harmonicvibro-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, andM. Fink, “Nonlinear shear wave interaction in soft solids,” J. Acoust.Soc. Am., vol. 122, no. 4, pp. 1917–1926, October 2007.71Chapter 3. Theoretical Limitations[20] A. Thitaikumar, T. A. Krouskop, B. S. Garra, and J. Ophir, “Visu-alization of bonding at an inclusion boundary using axial-shear strainelastography: 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. Tan-ter, and M. Fink, “Acoustoelasticity in soft solids: Assessment of thenonlinear 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, andJ. Ophir, “Breast tumor classification using axial shear strain elastog-raphy: a feasibility study,” Phys. Med. Biol., vol. 53, pp. 4809–4823,2008.[23] R. Lerner and K. Parker, “Sonoelasticity images, ultrasonic tissue char-acterization and echographic imaging,” in Proceedings of the 7th Euro-pean Communities Workshop, 1987, nijmegen, Netherlands.[24] R. Lerner, K. Parker, J. Holen, R. Gramiak, and R. Waag, “Sono-elasticity: Medical elasticity images derived from ultrasound signals inmechanically vibrated targets,” Acoustical Imaging, vol. 16, pp. 317–327, 1988.[25] R. Sinkus, J. Lorenzen, D. S. amd M. Lorenzen, M. Dargatz, andD. Holz, “High-resolution tensor mr elastography for breast tumor de-tection,” Physics Med. Biol., vol. 45, pp. 1649–1664, 2000.[26] R. Souchon, L. Soualmi, M. Bertrand, J.-Y. Chapelon, F. Kallel, andJ. Ophir, “Ultrasonic elastography using sector scan imaging and aradial compression,” Ultrasonics, vol. 40, pp. 867–871, 2002.[27] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulusimaging with 2-d transient elastography,” IEEE Trans. Ultrason. Fer-roelect. 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 elas-tography in anisotropic medium: Application to the measurement of72Chapter 3. Theoretical Limitationsslow 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 imag-ing: 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 andmotion 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-adaptivemotion tracking of ultrasound image sequences using a deformablemesh,” IEEE Trans. Med. Imaging, vol. 17, no. 6, pp. 945–956, 1998.[33] E. Konofagou, T. Varghese, J. Ophir, and S. Alam, “Power spectralstrain 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 ofnormal 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, “Di-rect 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 IEEEUltrasonics Symposium, 2002, pp. 1811–1820.[37] J. Bercoff, R. Sinkus, M. Tanter, and M. FinkJ, “3d ultrasound-baseddynamic and transient elastography: First in vitro results,” in IEEEUltrasonics Symposium, 2004, pp. 28–31.[38] K. Hoyt, F. Forsberg, and J. Ophir, “Investigation of parametric spec-tral estimation techniques for elasticity imaging,” Ultrasound Med.Biol., vol. 31, no. 8, pp. 1109–1121, 2005.[39] ——, “Comparison of shift estimation strategies in spectral elastogra-phy,” Ultrasonics, vol. 44, pp. 99–108, 2006.[40] ——, “Analysis of a hybrid spectral strain estimation technique in elas-tography,” Phys. Med. Biol. 51, vol. 51, pp. 197–209, 2006.73Chapter 3. Theoretical Limitations[41] J.-L. Gennisson, T. Deffieux, R. Sinkus, P. Annic, M. Pernot, F. Cud-eiro, G. Montaldo, M. Tanter, M. Fink, and J. Bercoff, “A 3d elastogra-phy system based on the concept of ultrasound-computed tomographyfor 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 noninvasivedetermination of material parameters from a knowledge of elastic dis-placements: 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 modu-lus 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 elas-ticity 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 linearperturbation 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 re-construction 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 inverseproblems in elasticity imaging using the adjoint method,” Inverse Prob-lems, vol. 19, pp. 297–313, 2003.[49] D. Fu, S. Levinson, S. Gracewski, and K. Parker, “Non-invasive quan-titative reconstruction of tissue elasticity using an iterative forward ap-proach,” Phys. Med. Biol. 45 (2000) 14951509, vol. 45, pp. 1495–1509,2000.[50] A. Manduca, T. Oliphant, M. Dresner, J. Mahowald, S. Kruse, E. Am-romin, J. Felmlee, J. Greenleaf, and R. Ehman, “Magnetic resonanceelastography: Non-invasive mapping of tissue elasticity,” Medical ImageAnalysis, vol. 5, pp. 237–2540, 2001.74Chapter 3. Theoretical Limitations[51] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complex-valued stiffness reconstruction by for magnetic resonance elastographyby algebraic inversion of the differential equation,” Magnetic Resonancein Medicine, vol. 45, pp. 299–310, 2001.[52] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelka-ram, and J. Culiolic, “Measurement of viscoelastic properties of ho-mogeneous soft solid using transient elastography: An inverse problemapproach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec2004.[53] B. Robert, R. Sinkus, B. Larrat, M. Tanter, and M. Fink, “A newrheological 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 withsurface 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, “Eval-uation 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 moduluselastography,” Ultrasound Med. Biol., vol. 31, no. 6, pp. 787–802, 2005.[57] M. M. Doyley, S. Srinivasan, E. Dimidenko, N. Soni, and J. Ophir, “En-hancing the performance of model-based elastography by incorporatingadditional a priori information in the modulus image reconstructionprocess,” 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 us-ing 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 andelastic properties ofsoft tissue using supersonic shear imaging,” in IEEEUltrasonics Symposium, 2003, pp. 925–928.[60] J. Fromageau, J.-L. Gennisson, C. Schmitt, R. L. Maurice, R. Mon-grain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel me-chanical properties with four ultrasound elastography methods and75Chapter 3. Theoretical Limitationscomparison with gold standard testing,” IEEE Trans. Ultrason. Fer-roelect. Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007.[61] A. Romano, J. Bucaro, P. Abraham, and S. Dey, “Inversion methodsfor the detection and localization of inclusions in structures utilizingdynamic surface displacements,” in Proceedings of SPIE Vol. 5503, vol.5503, 2004, pp. 367–374.[62] H. Kolsky, Stress Waves in Solids. New York: Dover Publications,Inc., 1963.[63] S. Catheline, F. Wu, and M. Fink, “A solution to diffraction biases insonoelasticity: The acoustic impulse technique,” J. Acoust. Soc. Am.,vol. 105, no. 5, May 1999.[64] A. Baghani, H. Eskandari, T. Salcudean, and R. Rohling, “Measure-ment of tissue elasticity using longitudinal waves,” in Proceedings ofthe 7th International conference on the Ultrasonic Measurement andImaging of Tissue Elasticity, October 2008, p. 103.[65] F. A. Duck, Physical Properties of Tissue: A Comprehensive ReferenceBook. Academic Press, 1990.[66] M. M. Burlew, E. L. Madsen, J. A. Zagzebski, R. A. Bahjavic, andS. 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 compilationof 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 medicaldiagnostics,” Ultrasound Med. Biol., vol. 24, no. 9, pp. 1419–1435, 1998.[69] M. Walz, J. Teubner, and M. Georgi, “Elasticity of benign and ma-lignant breast lesions,” in Imaging, Application and Results in Clinicaland General Practice, 1993, p. 56, eight International Congress on theUltrasonic 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 humanliver and correlation with pathology,” Ultrasound Med. Biol., vol. 28,no. 4, pp. 467–474, 2002.76Chapter 3. Theoretical Limitations[71] T. Krouskop, T. Wheeler, F. Kallel, B. Garra, and T. Hall, “Elasticmoduli of breast and prostate tissues under compression,” UltrasonicImaging, vol. 20, no. 4, pp. 260–274, October 1998.[72] M. Muller, J.-L. Gennisson, T. Deffieux, R. Sinkus, P. Annic, G. Mon-taldo, M. Tanter, and M. Fink, “Full 3d inversion of the viscoelasticitywave propagation problem for 3d ultrasound elastography in breastcancer diagnosis,” in IEEE Ultrasonics Symposium, 2007, pp. 672–675.[73] L. Huwart, F. Peeters, R. Sinkus, L. Annet, N. Salameh, L. C. terBeek, Y. Horsmans, and B. E. V. Beers, “Liver fibrosis: Non-invasiveassessment 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. Son-dermann, and M. Fink, “Imaging anisotropic and viscous propertiesof breast tissue by magnetic resonance-elastography,” Magnetic Reso-nance in Medicine, vol. 53, no. 2, pp. 372–387, 2005.[75] R. Sinkus, K. Siegmann, T. Xydeas, M. Tanter, C. Claussen, andM. Fink, “Mr elastography of breast lesions: Understanding thesolid/liquid duality can improve the specificity of contrast-enhanced mrmammography,” Magnetic Resonance in Medicine, vol. 58, pp. 1135–1144, 2007.[76] R. Baierlein, “Representing a vector field: Helmholtz’s theorem derivedfrom a fourier identity,” Am. J. Phys., vol. 63, pp. 180–182, 1995.[77] F. Rohrlich, “Causality, the Coulomb field, and Newton’s law of gravi-tation,” Am. J. Phys., vol. 70, no. 4, pp. 411–414, April 2002.[78] K. F. Graff, Wave Motion in Elastic Solids. Ely House, London:Oxford University Press, 1975.[79] L. Sandrin, D. Cassereau, and M. Fink, “The role of the coupling termin transient elastography,” J. Acoust. Soc. Am., vol. 115, no. 1, pp.73–83, January 2004.77Chapter 4The Influence of theBoundary Conditionson Longitudinal WavePropagationin a Viscoelastic Medium∗4.1 IntroductionThe propagation of an acoustic wave in soft tissue is governed by localmechanical properties. For tissue characterization, such as in the field ofelastography, the goal is to estimate the mechanical properties from mea-surements of the wave propagation. In such an inverse problem, the un-known mechanical parameters can be identified by analyzing the dynamicmotion of the medium. Of the two types of waves that can propagate ina viscoelastic environment, shear waves have received significantly more at-tention than longitudinal waves in the field of tissue characterization. Fora shear wave, particle displacements are in a direction perpendicular to thedirection of the wave propagation, while for a longitudinal wave, the parti-cle motions are in the same direction as the wave propagation. Shear wavesare usually produced inside tissue either by external transverse motion ofthe boundary [4,29], or through remotely induced acoustic radiation forceimpulse (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. Sal-cudean and R. Rohling, “The Influence of the Boundary Conditions on Longitudinal WavePropagation in a Viscoelastic Medium,” Phys. Med. Bio., vol. 54, no. 13, pp. 3997-4017,Jun. 200978Chapter 4. Influence of Boundary Conditionsultrasound or by magnetic resonance elastography (MRE) and the shearmodulus can be reconstructed from the velocity of the wave front [19].Longitudinal waves are assumed to have a very high propagation speedwhich makes them untraceable by ultrasound or other relatively low frame-rate data acquisition methods. Therefore, the effect of longitudinal wave isusually eliminated by applying the curl operator to the displacement field[25]. We have recently shown that it is possible to track the motion ofthe medium induced by longitudinal waves [1, 8]. Observations of wavepropagation in tissue-mimicking phantoms indicate that a longitudinal wavecan 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.Usingtheultrasoundimagingmodalityandalongitudinalexcitationscheme,the relaxation behavior of phantoms and breast tissue were observed andquantified in [12] and [13].In [22], a method for generating low-speed longitudinal waves was pro-posed. A low-frequency compressional excitation was applied by two rodsthat were placed on both sides of an ultrasound probe. As a result, a wavecould be seen in front of the probe, having the same speed as the shearwave. This phenomenon has been explained in [22] and [23] as the effect ofsuperposition of the shear wave fronts of the two exciters, superimposing inthe middle line and producing a longitudinally polarized shear wave throughmode conversion. Based on a Green’s function analysis in [21], this effectwas associated with the presence of a coupling term which gives rise to alongitudinal wave propagating at the shear wave speed.In this paper we analyze the effect of the boundary conditions on thepropagation speed of longitudinal waves in a finite medium. It will be shownhere by exploiting the three-dimensional (3D) wave equation that a longitu-dinal wave may travel at a very high speed such as ultrasound or at speedsas low as that of the shear wave. In fact, for nearly incompressible softmaterials with appropriate boundary conditions, a thin rod situation canbe approximated where a longitudinal wave travels at a speed of nearly 1.7times that of shear wave. The theoretical results will be verified by harmonicand transient 3D finite element simulations and the effect of viscosity willbe explored. A 3D linear dynamic FEM model is adopted as a referenceframework for the analysis and to validate the analytical results. Transientand harmonic simulations are performed. Experiments on a homogeneousincompressible soft material with different boundary conditions validate thelongitudinal wave speeds predicted by the analytical derivations and sim-ulations. As a result of this study, besides the shear wave, tracking thelongitudinal wave in ultrasound viscoelasticity imaging can be used as a79Chapter 4. Influence of Boundary Conditionsnew method to investigate the viscoelastic behavior of soft tissue.4.2 Governing Equations in a Linear ViscoelasticMediumIntheCartesiancoordinatesystem, the3Ddisplacementvector vectoru = [ux, uy, uz]Tand the strain tensor epsilon1 = [epsilon1xx, epsilon1yy, epsilon1zz, γxy, γxz, γyz]T have the followingrelationship:epsilon1 =∂x 0 00 ∂y 00 0 ∂z∂y ∂x 0∂z 0 ∂x0 ∂z ∂yvectoru . (4.1)In a linear isotropic viscoelastic medium, the stress tensor σ is relatedto the time-dependent strain tensor through a combination of the Lam´econstants λ and µ, and viscosity constants λprime and µprime:σ =λ +2µ λ λ 0 0 0λ λ +2µ λ 0 0 0λ λ λ +2µ 0 0 00 0 0 µ 0 00 0 0 0 µ 00 0 0 0 0 µepsilon1+λprime +2µprime λprime λprime 0 0 0λprime λprime +2µprime λprime 0 0 0λprime λprime λprime +2µprime 0 0 00 0 0 µprime 0 00 0 0 0 µprime 00 0 0 0 0 µprime∂tepsilon1 , (4.2)where ∂t is the time-derivative. Here, µ and µprime are often referred to asshear elasticity and shear viscosity, and λprime is known as the elongational orcompressional viscosity.The equation of motion for a body that is subjected to an external force80Chapter 4. Influence of Boundary Conditionsfield vectorf, is:ρ∂2t vectoru =∂xσxx +∂yσxy +∂zσxz∂xσxy +∂yσyy +∂zσyz∂xσxz +∂yσyz +∂zσzz+ vectorf , (4.3)where vectoruisthedisplacementvectorasafunctionoftimeandposition, andρisthe density. By inserting from (4.1) and (4.2) into (4.3), a hyperbolic partialdifferential equation (PDE) results which describes the wave propagation ina linear viscoelastic medium [16]:ρ∂2t vectoru = (λ+µ)∇∇·vectoru+µ∇2vectoru+(λprime +µprime)∇∇·(∂tvectoru)+µprime∇2(∂tvectoru)+ vectorf .(4.4)The gradient and divergence operators are denoted by (∇) and (∇·), re-spectively. ∇2 is the vector Laplacian operator. Note that in the Cartesiancoordinate system∇·vectoru = epsilon1xx+epsilon1yy+epsilon1zz, which is a measure of the volumetricchange in the medium and is called the dilatation.The Young’s modulus E is defined as:E = 2(1+ν)µ , (4.5)where,ν = λ2(λ+µ) . (4.6)ν is the Poisson’s ratio, which takes values between 0 and 0.5. A valueclose to 0.5 denotes incompressibility of the material, as is the case for mostsoft tissues. Different assumptions have been made in the literature to re-duce the number of viscosity parameters. For example, if the viscous effectsare assumed to follow the elastic effects in the material, then λprime/λ = µprime/µ.However, since there has been little justification for such similar behaviorof viscosity and elasticity, this assumption seems implausible [26]. Anotherassumption as made by [2] and [25] is that the compressional viscosity λprime isnegligible compared to the shear viscosity. To our knowledge, there is noexperimental validation for this assumption for incompressible soft tissue.Another approach is to translate the idea of elastic incompressibility into theviscosity domain by assuming that the bulk viscosity is negligible, so that thefluid incompressibility condition holds. When the divergence of flow is zero,the fluid is considered to be incompressible and the term associated withthe bulk viscosity in the Navier-Stokes equations will vanish [18]. In thiscase, the power dissipation in the medium is only a diviatoric dissipation81Chapter 4. Influence of Boundary Conditionsassociated with the distortion or shape-change, rather than dilatational dis-sipation that involves a change of volume [18]. The bulk viscosity is definedas κ = λprime +2µprime/3 which can be set to zero to obtain:λprime = −23µprime . (4.7)Bulk viscosity is known to be negligible compared to shear viscosity forrubber-like materials as well as some liquids and gases [16]. We proceed withthe 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 fordifferent sets of boundary conditions. The effective speed of the longitudinalwave propagation will then be extracted from the resulting representationof the wave equation.4.3 Wave Equation and the Propagation SpeedIn a purely elastic infinite medium, the wave equation (4.4) becomes:ρ∂2t vectoru = (λ+µ)∇∇·vectoru+µ∇2vectoru+ vectorf . (4.8)Thegeneralsolutionofequation(4.8)canbedecomposedintoadivergence-free and a curl-free component [16]. The divergence-free part tends to main-tain the volume and is referred to as the shear or s-wave which travels at aspeed of:cs =radicalbiggµρ . (4.9)The curl-free component has no rotation within infinitesimal blocks of mate-rial and is called the dilatational or p-wave which travels at a higher speed:cp =radicalBiggλ+2µρ . (4.10)While plane shear or dilatational waves in an infinite medium travel atthe designated speeds mentioned above, the influence from the boundaryconditions may give rise to waves that travel at other velocities. A fewexamples are presented in the following sections.82Chapter 4. Influence of Boundary Conditions4.3.1 Longitudinal Wave in a Finite MediumApproximating a Thin RodLet us first consider the case where the excitation is mainly compressional inthe x direction. If the exciter is large relative to the cross-section such thata plane-wave is produced in the medium, and the sides of the region havefree boundaries, all the stress components, except for σxx, are negligible.Note that for the viscoelastic case, the strain distribution will be affectedby the viscous stresses as well as the elastic stresses, however due to thecompressional excitation and the plane-wave assumption, σxx is dominantand the approximation should be valid. Using this condition in (4.2),bracketleftbig(λ+µ)+(λprime +µprime)∂tbracketrightbigσxx =parenleftbigµ+µprime∂tparenrightbigbracketleftbig(3λ+2µ)+(3λprime +2µ)∂tbracketrightbigepsilon1xx .(4.11)For the interior region that is not experiencing any internal force, (4.3) isreduced to ∂xσxx = ρ∂2t ux. This can be substituted into (4.11), togetherwith (4.1), to obtain the viscoelastic wave equation in the x-direction:ρbracketleftbig(λ+µ)+(λprime +µprime)∂tbracketrightbig∂2t ux = parenleftbigµ+µprime∂tparenrightbigbracketleftbig(3λ+2µ)+(3λprime +2µ)∂tbracketrightbig∂2xux .(4.12)This is a third order PDE with respect to the variable t. The case oflongitudinal wave propagation in a viscoelastic thin rod has been addressedmore than seven decades ago by Thompson [26]. The conditions described inthis section in which non-axial stress components become zero approximate athin rod situation. Therefore, following [26], by taking the Laplace transformin 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ωµprime)[(3λ+2µ)+jω(3λprime +2µprime)][(λ+µ)+jω(λprime +µprime)] = 0 . (4.13)There are three roots associated with this characteristic equation with re-spect to ω. However, considering that at relatively low frequencies (λ+µ) greatermuch(λprime+µprime)ω for incompressible materials, the root that corresponds to the de-nominator is very large and may be ignored [26]. Also, with the assumptionof zero bulk viscosity, 3λprime +2µprime = 0. Therefore, equation (4.13) reduces to:−ω2ρ+s2(µ+jωµprime)(3λ+2µ)(λ+µ) = 0 . (4.14)83Chapter 4. Influence of Boundary ConditionsGiven the incompressibility of the material and the definition of Young’smodulus, this can be written as −ω2ρ + s2(E + jωη) = 0, which results inthe following PDE:ρ∂2t ux = E∂2xux +η∂2x∂tux . (4.15)η is called the coefficient of normal viscosity, where,η ≈ 3µprime . (4.16)The PDE in equation (4.15) can be discretized using the well-known dis-tributed system of several mass-spring-damper (MSD) elements which areconnected to each other in series [27]. Based on (4.16), in order for the 1Dmodel to follow the behavior of the 3D model in its best sense, the equiv-alent viscosity of the 1D model should be three times the shear viscosity.This will be taken into account in the following simulations where behaviorsof different models are compared. Also, given that the Young’s modulusin nearly incompressible materials is approximately three times the shearmodulus, the ratio between the observed viscous and elastic effects, enti-tled 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 themedium, propagating at a much lower speed than cp. In the purely elasticcase, without the influence of viscosity, the speed is:c0 =radicalBiggEρ . (4.17)This value, being approximately 1.7 times the shear wave speed for in-compressible materials, denotes that a low-speed longitudinal wave can beinduced in the medium by applying certain boundary conditions. Notethat if the viscosity is also taken into account, the effective modulus from(4.15) will be radicalbigE2 +ω2η2 instead of E, which results in a speed equal toc ≈ c0(1 + ω2τ2/4), where τ = η/E is the relaxation-time. However, withthe small viscosity and relaxation-time in soft tissue, this correction termis negligible at low frequencies and equation (4.17) can be used to calculatethe wave speed.4.3.2 Longitudinal Wave Induced by a Point Source ExciterSuppose a point source excitation is applied axially to a finite sample ofhomogeneous material in which radial symmetry can be assumed. In this84Chapter 4. Influence of Boundary Conditionsuruè uxrèxfxFigure 4.1: Displacement vector in the cylindrical coordinate system.case, it is more convenient to analyze the problem in cylindrical coordi-nates as shown in Fig. 4.1. The stress-strain relationship in the cylin-drical coordinate system is the same as (4.2). Similar to the analysis inSection 4.2, in cylindrical coordinates where vectoru = [ux, ur, uθ]T and epsilon1 =[epsilon1xx, epsilon1rr, epsilon1θθ, γxr, γxθ, γrθ]T, the strain tensor is defined as follows:epsilon1 =∂x 0 00 ∂r 00 1r 1r ∂θ∂r ∂x 01r ∂θ 0 ∂x0 1r ∂θ (∂r −1r)vectoru . (4.18)The equation of motion in the x-direction in this case becomes:ρ∂2t ux = ∂xσxx +∂rσxr + 1rσxr + 1r ∂θσxθ +fx . (4.19)Note that if all of the shear stress components are zero, the resulting waveequation would be that of the p-wave in a viscoelastic medium.For the interior, where fx = 0, (4.19) can be expanded as:ρ∂2t ux = [(λ+2µ)+(λprime +2µprime)∂t]∂x(epsilon1xx +epsilon1rr +epsilon1θθ)−2(µ+µprime∂t)∂x(epsilon1rr +epsilon1θθ)+(µ+µprime∂t)∂repsilon1xr + 1r(µ+µprime∂t)epsilon1xr + 1r(µ+µprime∂t)∂repsilon1xθ .(4.20)85Chapter 4. Influence of Boundary ConditionsIt can be shown that equation (4.20) can be simplified to:ρ∂2t ux = [(λ+2µ)+(λprime +2µprime)∂t]∂x(epsilon1xx +epsilon1rr +epsilon1θθ)−2(µ+µprime∂t)1r [∂r(r¯ωθ)−∂θ(¯ωr)] , (4.21)where ¯ωr and ¯ωθ are rotations about the r and θ axes:2¯ωθ = ∂xur −∂rux , and, 2¯ωr = 1r ∂θux −∂xuθ . (4.22)To obtain a relationship between the radial and angular strains on theaxis of excitation, l’Hˆopital’s rule can be utilized at r → 0. According tothis theorem, if two functions are differentiable in a neighborhood of a givenpoint while both of them converge to zero, the indeterminate limit of theirquotients would be equal to the quotient of their derivatives at that point.By applying this rule to the angular strain, one gets:epsilon1θθ|r=0 = limr→0 1rur = limr→0∂rur = epsilon1rr|r=0 . (4.23)Also, due to radial symmetry in this case, σθθ = 0, which according to(4.2), leads to:bracketleftbig(λ+2µ)+(λprime +2µprime)∂tbracketrightbigepsilon1θθ +(λ+λprime∂t)(epsilon1xx +epsilon1rr) = 0 . (4.24)Inserting from (4.23) into (4.24), at the axis of excitation:epsilon1rr = epsilon1θθ = −(λ+λprime∂t)2[(λ+µ)+(λprime +µprime)∂t] epsilon1xx . (4.25)Inserting (4.25) into (4.21), noting that ¯ωr is zero on the excitation axisand assuming that ¯ωθ and ∂rˆωθ are also zero at r → 0, the following results:ρbracketleftbig(λ+µ)+(λprime +µprime)∂tbracketrightbig∂2t ux = parenleftbigµ+µprime∂tparenrightbigbracketleftbig(λ+2µ)+(λprime +2µ)∂tbracketrightbig∂2xux .(4.26)With the same approach as in section 4.3.1 and given the large value ofλ compared to µ,µprime and λprime in incompressible materials, the wave equationsimplifies to:ρ∂2t ux = µ∂2xux +µprime∂2x∂tux . (4.27)Therefore, given radial symmetry in the medium, a point-source exciter86Chapter 4. Influence of Boundary Conditionswill generate a longitudinal wave that travels at the shear speed:cprime0 =radicalbiggµρ , (4.28)The assumption that ¯ωθ and ∂rˆωθ vanish at r → 0 is verified in Section 4.5 byobserving the same wave as predicted above in simulations and experiments.Also, note that using (4.25), the axial stress will be:σxx = [(λ+2µ)+(λprime +2µprime)∂t]epsilon1xx +(λ+λprime∂t)(epsilon1rr +epsilon1θθ)= (µ+µprime∂t)[(3λ+2µ)+(3λprime+2µprime)∂t][(λ+µ)+(λprime+µprime)∂t] epsilon1xx≈ Eepsilon1xx +η∂tepsilon1xx , (4.29)where η is as in equation (4.16). In other words, the apparent elastic con-stant in the axial direction is the Young’s modulus rather than the shearmodulus, however the velocity of the wave is determined by the shear mod-ulus according to (4.28).In the following section, it will be shown that for the case of a homo-geneous block of material, the displacement spectra and particularly theresonance frequency of the system can be used to estimate the wave speedin the medium. The effect of viscosity on the resonance frequency and thespeed will be formulated. The analysis in the following section will justify amethod that will be used in Section 4.5 in order to estimate the wave speedfrom the displacement spectrum.4.4 Determining the Wave Speed from theSpectra of the DisplacementsIn this section, the effect of the viscosity on the wave equation, and in partic-ular on the displacement spectrum, will be analyzed. It will be shown thatthe wave speed can be estimated from the spectrum of the displacements ina viscoelastic medium. The result from this section will be used to estimatethe wave speed from the simulation results in Section 4.5, where the effect ofthe boundary conditions on the wave speed is explored. As shown in [1], ifthe medium is bounded from one end, the solution to (4.15) in the frequencydomain is in the following form:ˆu(x,ω) = A0 [cos(βx)sinh(αx)+jsin(βx)cosh(αx)] , (4.30)87Chapter 4. Influence of Boundary Conditionswhere ˆu is the Fourier transform of ux, A0 is the amplitude of the vibration atfrequency ω and α and β can be calculated by taking the Laplace transformof (4.15) in the x domain and the Fourier transform in the t domain [1]:(E +jωη)s2 +ω2ρ = 0 . (4.31)Definingrelaxation-time byτ, insofttissue, ωτ = ωη/E lessmuch 1forω < 100Hz,the solution to (4.31) can be written as:s = ±(α+jβ) = ± jωc0√1+jωτ ≈±jωc0(1− jωτ2 ) , (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 frequen-cies, it can be shown through simple trigonometric and hyperbolic transfor-mations, 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 speedc0 can be obtained. Two points on this functions are of particular interestto us:(i) The frequency f1 at which the displacement at x = L/2 is largercompared to other locations.(ii) The frequency f2 at which the displacement at x = L/2 has itsmaximum 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 frequencythat maximizes the displacement at the middle point or the frequency forwhich the middle point has the highest displacement magnitude comparedto 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 highviscosity medium. A Young’s modulus of 10kPa and density of 1000kg/m3were assumed, thus c0 = 3.16m/s and f0 = c0/(2L) = 26.4Hz. First, with88Chapter 4. Influence of Boundary Conditionsa viscosity of 5Pas, α and β were calculated from (4.32). The displacementmagnitude 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 f1and f2 respectively. It can be seen that if the viscosity is very low, the twocases are almost the same and both f1 and f2 occur very close to f0. Inthe second case, a viscosity of 30Pas was assumed and the displacementmagnitudes were calculated. The magitude and contour plots are depictedin 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 thefollowing, 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 approximationproduces will be quantified. The more accurate method will be adopted forwave 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,12ω1τsinh(αL)+sin(βL) = 0 . (4.36)Basedontheassumptionthatωτ lessmuch 1andthatωL lessmuch c0, wehavesinh(α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 ≈ pi. Thus, sin(βL) ≈ pi −βL. From (4.36)and (4.33):ω31τ2L4c0 +pi −ω1Lc0 = 0 . (4.37)If ω1 = 2pif1, the wave speed can be estimated as follows:c0 = 2f1Lparenleftbig1−pi2f21τ2parenrightbig . (4.38)Apparently in this case, if one uses 2f1L instead of (4.38), the wave speedwill be over-estimated by a small fraction.89Chapter 4. Influence of Boundary ConditionsFrequency (Hz)Depth (cm)1020304050600 1 2 3 4 5 61020304050(a)Frequency (Hz)Depth (cm)2040606543210f1, f2(b)Frequency (Hz)Depth (cm)1020304050600 1 2 3 4 5 60.511.5(c)Frequency (Hz)Depth (cm)2040606543210f1f2(d)Figure 4.2: The magnitude of the displacement in a 1D model of the waveequation with the purely elastic wave speed of 3.16m/s, length of 6cm, andYoung’s modulus of 10kPa. (a) and (b) are the magnitude and contour plotsfor the case in which the viscosity is 5Pas. The frequencies denoted by case(i) and case (ii) are almost equal. (c) and (d) show the the plots for the casethat the viscosity of the medium is 30Pas. In this case, f1 over-estimatesthe 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 itto zero for x = L/2, the following will be obtained:(ω2τ sinh(αL)+sin(βL))parenleftbigsinh2(αL)+sin2(βL)parenrightbig =2(ω2τ sinh(2αL)+sin(2βL))parenleftbiggsinh2(αL2 )+sin2(βL2 )parenrightbigg. (4.39)90Chapter 4. Influence of Boundary ConditionsUsing the characteristics of the trigonometric and hyperbolic functions, itcan be shown that:sinh2(αL)+sin2(βL) = (cosh(αL)+cos(βL))(cosh(αL)−cos(βL)) ,(4.40)and,sinh2(αL2 )+sin2(βL2 ) = cosh(αL)−cos(βL)2 . (4.41)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 becalculated as:c0 = 2f2Lparenleftbig1+pi2f22τ2parenrightbig . (4.44)where ω2 = 2pif2. If a small amount of under-estimation is tolerable, thewave speed can be approximated as 2f2L.By comparing (4.38) and (4.44), it can be seen that f2 < f1, therefore thesecond-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 thewave speed as 2f2L rather than 2f1L, if the relaxation-time of the mediumis unknown. Also, note that although equation (4.15) has been used here toanalyze the accuracy of the two cases above, the same argument is valid forother forms of the wave equation, for example when shear modulus replacesYoung’s modulus as in (4.27).4.5 Results4.5.1 Dynamic Finite Element AnalysisNumerical analysis of the three-dimensional (3D) wave equation is possiblethrough discretization of equation (4.4) using the Ritz-Galerkin method [6].The resulting finite element model can be used to explain the static and dy-namic deformations in the medium in a discretized platform. An equationcan be obtained to describe the relationship between the nodal displace-91Chapter 4. Influence of Boundary ConditionsXZY4cm6cm4cm[0,0,0]TExcitationFigure 4.3: The geometry of the 3D region used for FEM simulations.ments, boundary conditions and local material properties:Ku(t)+B˙u(t)+M¨u(t) = f(t) , (4.45)where K, B and M are the stiffness, viscosity and mass matrices of theelements, u(t) is the nodal displacement vector and f(t) is a vector contain-ing the boundary conditions. This equation can be written in the frequencydomain as follows:(K +jωB−ω2M)ˆu = ˆf . (4.46)A method to generate an accurate viscosity matrix based on a linearviscoelastic model of the continuum has been proposed in [9]. An implicitsolution for the transient case of (4.45) with step-size ∆t can be obtainedas [3]:¯Ku(t+∆t) = ¯f(t) ,¯K = K + B∆t +M∆t2 ,¯f(t) = f(t)+ 2M∆t2u(t)+(B2∆t −M∆t2)u(t−∆t) . (4.47)A finite element simulation can be used to verify the theory presentedin Section 4.3. In particular, numerical tests can facilitate analysis of thepropagation of the fast and slow longitudinal waves in soft materials. Forthis purpose, a region of size 6×4×4(cm3) was meshed by eight-node brickelements, having a total of 25×7×7 nodes. The bottom surface at x = 0was confined from axial movement in the x direction and the excitation92Chapter 4. Influence of Boundary Conditionswas applied axially at the top surface, as shown in Fig. 4.3. Although thederivations in section 4.3 have been based on either radial symmetry orthin rod assumptions, the Cartesian-based rectangular FEM model utilizedhere does not satisfy those conditions. However, the geometry used herecan provide additional insight into the longitudinal wave propagation in anon-ideal situation.In this manuscript, the middle line refers to the line in the phantomwith {y,z} = [2,2]cm, the middle point or middle node refers to the pointor node on top of the phantom with {x,y,z} = [6,2,2]cm, and the centerpoint is the point at {x,y,z} = [3,2,2]cm.4.5.2 Simulation of the Wave PropagationThe 3D region in Fig. 4.3 was assumed to consist of a homogeneous materialwith a Young’s modulus of 10kPa, a viscosity of 10Pas and a density of1000kg/m3. A Poisson’s ratio of ν = 0.495 is often considered for finiteelement analysis of soft tissues [11]. With the given Young’s modulus andPoisson’s ratio, it can be shown that λ = 334.1kPa and µ = 3.3kPa.Three settings for the boundary conditions were simulated. The longi-tudinal wave patterns were compared to the theoretical results in Section4.3 for the following cases, when:(a) The excitation was applied to all of the nodes at the top surface ofthe phantom. The four vertical sides of the phantom were free, thusbulging was allowed.(b) The excitation was applied to all of the nodes at the top surface ofthe phantom. The four vertical sides of the phantom were confinedwith slip boundaries such that bulging was not allowed.(c) The excitation was applied to a single node at the center of the topsurface, i.e. the middle node. The four vertical sides of the phantomwere free, thus bulging was allowed.In all these cases, a slip boundary condition was imposed on the bottomof 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 possibledue to the boundary conditions. Therefore, a dilatational wave will begenerated at a speed determined by (4.10). Finally, in case (c), a pointsource compressional exciter is modeled within the resolution of the finiteelement mesh. As shown in section 4.3.2, the speed of the longitudinalcomponent of the generated wave in this case is approximated by (4.28).These analytical derivations will be tested in the following simulations.93Chapter 4. Influence of Boundary Conditions10110210−310−210−1100101Frequency (Hz)|H(f)|TF Magnitude at the Center PointFree sidesBounded sidesPoint source(a)101102−6−4−20Frequency (Hz)∠ H(f) [rad]TF Phase at the Center PointFree sidesBounded sidesPoint source(b)Figure 4.4: The magnitude (a) and the phase (b) of the transfer functionbetween 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 phantomare free to move and the dashed lines are for the case where the sides ofthe phantom are confined from bulging. The gray dash-dot line shows thetransfer function for a one-node excitation at the middle point.For each case, a frequency sweep analysis was performed by solving equa-tion (4.46) over a range of frequencies. The transfer function between thedisplacements at the center point and the middle point was then measuredby dividing the displacements spectra at the two locations. As explained insection 4.4, by determining the frequency at which the peak of the magni-tude of this transfer function occurs, the wave speed can be estimated giventhe small relaxation-time of the medium.The transfer function magnitudes for the three cases above are shown inFig. 4.4(a) and the phases are plotted in Fig. 4.4(b). The peak of the threetransfer functions occur at 24.4Hz, 153Hz and 14.6Hz for cases (a), (b)and (c), respectively. These frequencies can be multiplied by a wavelengthof 12cm, which is twice the axial extent of the phantom, to obtain thelongitudinal wave speed. Therefore, the estimated wave speed for thesecases would be respectively, 2.93m/s, 18.36m/s and 1.75m/s. The actualvalues 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 mediumwould change. In case (a), the coefficient of normal viscosity can be cal-culated from (4.16) to be 30Pas, which should be divided by the Young’smodulus to get a relaxation-time of τa = 1ms. In case (b), the effective vis-cous modulus would be λprime+2µprime, which is a counterpart of the effective elastic94Chapter 4. Influence of Boundary Conditions10110210−1100101Frequency (Hz)|H(f)|Effect of Viscosity on the Magnitude of TF1 Pas5 Pas10 Pas20 Pas100 PasFigure 4.5: Effect of the viscosity on the magnitude of the transfer functionbetween the center point and the middle point. The excitation was appliedto the entire top surface while the sides were free to bulge.modulus λ+2µ. Given a bulk viscosity of zero, λprime+2µprime = 4µprime/3 = 13.3Pas.Thus, the relaxation-time in this case would be τb = 39µs. Note that in real-ity where λ is in the order of 109, the viscosity for p-wave propagation wouldbe even less significant. Lastly, for case (c), the viscosity would be equal toµprime according to (4.28), which can be divided by the shear modulus to geta 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 thetransfer function magnitude. Higher viscosity results in higher dissipationand less power transfer in the medium. Figure 4.5 shows the effect of chang-ing the viscosity on the magnitude of the transfer function for the samemodel as explained above. The elastic wave was generated in this case bycompressing the top of the phantom and allowing the sides to move freelyas in case (a). Note that the location of the peak and thus the speed ofpropagation do not change significantly, while the attenuation depends onthe damping value.To show the transient behavior of the model, a step displacement with1mm amplitude was applied to the top of the simulated phantom. Figure 4.6shows how the step propagates along the middle line of the phantom. Theprevious three cases were explored here as well. First, an elastic wave wasgenerated in a medium with free sides by applying the excitation to the topsurface. Based on the peak of the displacement waveform in Fig. 4.6(a), it95Chapter 4. Influence of Boundary Conditionstakes 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), theresonance of the phantom in response to the step is at 24.4Hz which complieswith the previous result. The magnified version of this figure is shown inFig. 4.6(b) for a few nodes on the middle line. This figure clearly illustratesthe way that information travels in the medium. It can be seen that a fasttraveling wave propagates in the phantom, which reaches the bottom after2.8ms. This corresponds to a speed of 21.4m/s which is comparable tothe p-wave speed in the medium, i.e. 18.38m/s. It should be noted thatdetermining the temporal location of the peak only gives an estimate ofthe wave-front rather than rendering it accurately. For the second test asin case (b), a compressional step excitation was applied to the top surfacewhile 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 themedium, a magnified version in Fig. 4.6(d) only shows the presence of onlyone wavefront at a speed of 18.2m/s. Finally, the middle node on the top ofthe phantom was excited by a step function. The middle line displacementsare depicted in Fig. 4.6(e). While the elastic wave speed is estimated at2.1m/s, a p-wave can be identified in the magnified version in Fig. 4.6(f) witha speed of 21.4m/s. Note that in all cases, the resonances that are producedin the phantom are consistent with the peaks of the transfer functions inFig 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 transferfunction for the center node, before and after viscosity compensation areshown. Also, the wave speed that is recovered from the temporal peak ofthe wave-front are reported for the three cases of boundary conditions.4.5.3 ExperimentsA soft homogeneous PVC (polyvinyl chloride) phantom was used to ex-perimentally validate the speeds that were derived in Section 4.3. Theshape of the phantom was as depicted in Fig. 4.3, but the dimensions were4 × 3 × 3(cm3) where the axial extent was 4cm. Independent rheometrymeasurements resulted in a Young’s modulus of 17kPa and a density of1000kg/m3 as reported in [1]. This will result in an elastic longitudinalwave speed of 4.1m/s and a shear wave speed of 2.4m/s. An ultrasoundlinear probe and a Sonix RP machine (Ultrasonix Corp., Richmond, BC,Canada) were used to capture RF data from a single A-line in the middleof the phantom at 5000 frames per second in the Doppler mode. The sides96Chapter 4. Influence of Boundary Conditionsof the phantom were free to bulge, the transducer was placed at the bottomand a step excitation was applied from the top. The average strain in thephantom due to the step was approximately 0.1%. The axial displacementswere estimated from the RF data using the peak search algorithm [7].In the first experiment, the exciter was a plate entirely covering the topsurface of the block. The displacements of different nodes along the acqui-sition line are shown in Fig. 4.7(a) where every displacement waveform isshifted up to the initial position of the block to which it corresponds. Byfollowing the peak of the wavefront in time, an estimate of the speed canbe obtained. The peak location of the wavefront is plotted versus time inFig. 4.7(b) for the first half of the phantom from the top. Fitting a line tothe data reveals a speed of 4.0m/s for the wave travelling longitudinally inthis case. One can also estimate the speed from the resonance frequency ofthe phantom. As seen in Fig. 4.7(a) the period of the vibration is 19.0mswhich denotes a resonance frequency of 52.6Hz. Given that this frequencycorresponds to a wavelength twice as large as the axial length of the phan-tom, 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 usinga small screw with a cross-section of 2mm2 to vibrate the phantom. Thedisplacement waveforms are shown in Fig. 4.7(c), where the initial frontTheory Estimation (m/s)Excitation (m/s) Spectrum withSpectrum viscosity Transientpeak compensation step(a) Flat excitation 3.16 2.93 2.95 2.9Free sides(b) Flat excitation 18.38 18.36 18.37 18.2Bounded(c) Point source 1.83 1.75 1.76 2.1Table 4.1: The theoretical and estimated wave speed for three differentconditions. The estimated wave speed are based on three method: first, thefrequency that the peak of the spectrum happens was used to reconstructthe wave speed; second, the previous value was enhanced to account for theeffect of viscosity; third, the wave-front due to a transient step was analyzedto recover the speed.97Chapter 4. Influence of Boundary ConditionsExcitation Theory Estimation (m/s)(m/s) Wavefront peak ResonanceFlat plate 4.1 4.0 4.2Point source 2.4 2.3 2.2Table 4.2: The theoretical and estimated wave speed for an experiment witha soft PVC phantom. A flat and a point-source excitation were applied toexplore the effect of boundary conditions on the speed of longitudinal waves.is delineated by a strain line passing through the peaks. A least-squaresline has been fitted to the peak locations for the onset of the wave in thephantom as shown in Fig. 4.7(d). The slope of this line indicates that alongitudinal wave travels at a speed of 2.3m/s in this case. Also, from thenatural vibration of the phantom, a period of 35.8ms can be infered whichcorresponds to a resonance frequency of 27.9Hz. Using the same calculationas for the previous case, the speed can be estimated at c = 2.2m/s. Theresults of these two experiments are summarized in Table 4.2.4.6 DiscussionThe speed of propagation of longitudinal waves is highly dependent on theboundary conditions as well as the material properties. In nearly incom-pressible materials, when the boundary conditions force the wave to propa-gate only through a change of volume, a high speed (cp) p-wave occurs. Thep-wave, traveling as fast as ultrasound in the medium, cannot be trackedby means of ultrasound, hence compressional transient elastography has re-ceived very limited attention in the literature. However, it is shown herethat with proper boundary conditions and excitation, it is feasible to inducea longitudinally propagating wave in the medium that has a much lowerspeed than the p-wave. Specifically, it is shown that if a block of materialcan freely preserve its volume, having a large uniform exciter on the sur-face can approximate the semi-infinite thin rod situation in which an elasticwave travels at only a few meters per second (c0). On the other hand, ifthe large exciter is replaced by a point-source indenter, the speed can benearly reduced to the shear wave speed in the medium (cprime0). It was pre-viously shown in [21], by using a Green’s function approach, that besides98Chapter 4. Influence of Boundary Conditionsthe s- and p-components of the wave, a coupling term is also generated ina semi-infinite medium which travels longitudinally in the near field at thespeed of the shear wave. However, in this paper we elaborated that thisspeed depends on the dimensions of the source of excitation relative to themedium dimensions. One can expect that a mid-size exciter can produceother velocities between cprime0 and c0 while tuning the boundary conditions andthe 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 changewould be inevitable and therefore, it is only the p-wave that conveys theinformation. Due to near incompressibility of the material, the p-wave gen-erated in the medium will be much faster than the longitudinal elastic wavewhen the sides are free.Finite element analysis was performed to validate the theoretical find-ings and explore the effect of the boundary conditions. To achieve stablefinite element solutions, the Poisson’s ratio is chosen smaller than the actualvalue. Although a closer value to 0.5 is more realistic for incompressible softtissues, it causes locking of the elements and instability in solving equation(4.45). Locking or hourgalssing of the elements make them behave in a stifferfashion. The problem of locking can be reduced by numerical integration in-stead of closed-form integration [6]; however, this may introduce inaccuracyin the final solution. A remedy to this trade-off is to assume near incom-pressibility, thus ν = 0.495. As a result, locking would not happen, howeverthe speed of the p-wave would be smaller in the simulations compared tothe actual speed of ultrasound in incompressible materials where Poison’sratio is closer to 0.5. A closer value of Poisson’s ratio to 0.5 results in highervalues of λ and cp. For example, having E = 10kPa and ν = 0.5−7×10−7gives cp = 1540m/s.The peak of the transfer functions has been used to calculate the wavespeed. At the resonance frequency of a purely elastic phantom, all the pointsvibrate at the same phase. The first mode corresponds to the frequency thata full wavelength is trapped in the finite medium, therefore the wavelengthof 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 highestamplitude of vibration which corresponds to a peak in the correspondingtransfer function magnitude as shown in Fig. 4.4(a). The phase also dropsatthenaturalfrequencyasseeninFig.4.4(b). Thesignificantphasedecreasein the simulations is a result of using a linear model, while real tissue oftenacts in nonlinear ways. As an example, human tissue and gelatin-basedphantoms are known to have a power-law viscous behavior, due to which99Chapter 4. Influence of Boundary Conditionsviscosity and relaxation-time drop as the frequency increases [8,14,15]. Thisnonlinearity causes the phase to have a smaller change at high frequencies.However, a linear FEM model is unable to account for such effects. Asshown in Fig. 4.5, having a lower viscosity gives rise to various resonancesin the system at high frequencies. Therefore, generally one may see severalpeaks in the transfer function, given a wide-band excitation.It is known that when a uniform quasi-static compression is appliedto a block of material from one side while the other side is bounded, theaxial displacements attenuate almost linearly. However, in the case of apoint-source exciter, the intensity falls off similarly to an ultrasound beamas the square of the distance from the source [5]. Therefore, the quasi-static displacement amplitude for the case of a point-source is expectedto be less than that of a large exciter. This explains the smaller transferfunction magnitude for a point-source compared to the other two cases inFig. 4.4(a) and the non-uniform spacing of the nodal displacement curves inFig. 4.6(e) while Figs. 4.6(a) and 4.6(c) exhibit uniform variation betweenthe displacements at the tail of the curves.It has been shown that a longitudinal wave may travel at a slow speedbetween c0 and cprime0 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 speedelastic wave is not the only wave present in the medium. Upon exertingan axial compression to a finite block, a p-wave is first generated in themedium which travels at the high speed of cp. At this instant, the wavehas not reached the boundaries and therefore, the boundary conditions andgeometry are not known. Once such information is collected by the p-wave,depending on the boundary conditions, a slower wave starts to propagatethat is influenced by the material properties as well as the boundary condi-tions. However, if the size of the block is infinite compared to the exciter -such as in ultrasound or acoustic waves - bulging and rotation of infinitesi-mal blocks would be infeasible. In this case, the irrotational wave propagatessolely through volumetric change of the medium at a speed close to cp.As noted before, when dealing with living tissue or tissue-mimickingmaterials where Poisson’s ratio is closer to 0.5, the speed of the p-waveis cp ≈ 1540m/s, while the speed of the slower longitudinal wave is notvery 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 inFig. 4.6(c). Note that in this case, the resonance of the p-wave lasts for alonger 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 themedium.100Chapter 4. Influence of Boundary ConditionsThe analytical and simulation findings were also validated experimen-tally on an incompressible viscoelastic phantom with known elastic proper-ties in Fig. 4.7. Only the onset of the wave front in the medium has beenused to estimate the wave speed. The reason is that the excitation was ap-plied with a voice-coil system which could not produce a perfect step withan infinitely small rise-time due to its open-loop control. As a result, bythe time the motion at the exciter has reached its peak, a wave had alreadystarted to propagate in the medium. This wave may reflect from the otherside of the block and interfere with the forward-traveling wave front some-where in the middle of the phantom depending on the speed of the wave.Due to this interference, the phase velocity changes beyond 20mm deep inthe phantom in Fig. 4.7(a) and beyond 15mm in Fig. 4.7(c). How far thepeak of the wavefront travels before this interference occurs depends on therise-time of the excitation, wave speed and the length of the medium.4.7 ConclusionIt is feasible to use longitudinal waves in order to identify the mechani-cal properties of soft tissues. Depending on the boundary conditions, alow speed longitudinal wave can be generated in the medium which canbe tracked by means of standard pulse-echo ultrasound. The speed of thiswave is determined by the material properties, exciter size and other bound-ary conditions. The longitudinal wave speed for a viscoelastic homogeneousmaterial can be found from the resonance frequency of the phantom, peakof the transfer function between a pair of locations, or phase speed of thetransient wavefront. The analytical findings were validated numerically us-ing 3D viscoelastic FEM models and experimentally on a tissue-mimickingphantom.101Chapter 4. Influence of Boundary Conditions02040608010000.20.40.60.81Time (ms)Displacement (mm)Unbounded Elastic Excitation(a)0510152025Time (ms)Displacement LogarithmUnbounded Elastic ExcitationExcitationSecond nodeCenter nodeLast nodeBottom node(b)02040608010000.20.40.60.81Time (ms)Displacement (mm)Bounded Medium(c)051015202500.20.40.60.81Time (ms)Displacement (mm)Bounded Medium(d)02040608010000.20.40.60.81Time (ms)Displacement (mm)Point−source Excitation(e)0510152025Time (ms)Point−source ExcitationDisplacement Logaorithm(f)Figure 4.6: Displacements at the middle line for a step excitation. The leftcolumn has a larger time-span to show the low-speed wave, where each linerepresents the displacement of one node. The high-speed wave can be seenin the short time-span of the right column. In (a) and (b) the excitation wasapplied to the top surface when the four sides were free. (c) and (d) showthe response for a single-node exciter at the middle node. In (e) and (f) themedium was bounded and could not bulge. Note that the displacement axesin (b) and (d) are in a logarithmic scale.102Chapter 4. Influence of Boundary Conditions05010015020001020304050Time (ms)Initial Position (mm)19.0 ms(a)1314151617184035302520Time (ms)Wavefront Peak Poistion (mm)|slope| = 4.0 m/sWave travel from the plate Fitted linear model(b)050100150200010203040Time (ms)Initial Position (mm)35.8 ms(c)1112131415403530Time (ms)Wavefront Peak Position (mm)|slope| = 2.3 m/sWave travel from the point−source Fitted linear model(d)Figure 4.7: Wave propagation in a soft PVC phantom with E = 17kPaand an axial length of 4cm under two different boundary conditions whena step excitation is applied. (a) and (c) show the displacements of tissueat different depths for a flat and a point-source excitation respectively. Thelocation of the wavefront peaks are shown versus time for these two boundaryconditions in (b) and (d), where a line has been fitted to the data to estimatethe speed. The resonance frequency of the phantom is also a function of theelastic modulus, geometry and boundary conditions.103References[1] Baghani, A., Eskandari, H., Salcudean, S. E. and Rohling, R.: 2009,Measurement of Viscoelastic Properties of Tissue Mimicking MaterialUsing Longitudinal Wave Excitation, Submitted to the IEEE Transac-tions on Ultrasonics, Ferroelectrics and Frequency Control.[2] Bercoff, J., Tanter, M. and Fink, M.: 2004, Supersonic shear imag-ing: a new technique for soft tissue elasticity mapping, IEEE TransUltrasonics, 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., Abouelka-ram, S. and Culioli, J.: 2004, Measurement of viscoelastic proper-ties of homogeneous soft solid using transient elastography: An in-verse problem approach, Journal of the Acoustical Society of America116(6), 3734–3741.[5] Cobbold, R. S. C.: 2007, Foundations of Biomedical Ultrasound, OxfordUniversity Press, NY.[6] Cook, R. D., Malkus, D. S. and Palesha, M. E.: 1989, Concepts andApplications of Finite Element Analysis, 3rd edn, John Wiley & Sons,Inc., New York.[7] Eskandari, H., Salcudean, S. E. and Rohling, R.: 2007, Tissue StrainImaging Using a Wavelet Transform-Based Peak Search Algorithm,IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Con-trol 54(6), 1118-1130.[8] Eskandari, H., Salcudean, S. E. and Rohling, R.: 2008a, Viscoelasticparameter estimation based on spectral analysis, IEEE Transactions onUltrasonics, Ferroelectrics and Frequency Control 55(7), 1611–1625.104Chapter 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 elementmodels, Physics in Medicine and Biology 53(22), 6569-6590.[10] Fatemi, M. and Greenleaf, J. F.: 1999, Vibro-acoustography: An imag-ing modality based on ultrasound-stimulated acoustic emission, Pro-ceedings 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 it-erative forward approach, Physics in Medicine and Biology45(6), 1495–509.[12] Insana, M. F., Pellot-Barakat, C., Sridhar, M. and Lindfors, K. K.:2004, Viscoelastic imaging of breast tumor microenvironment with ul-trasound, Journal of Mammary Gland Biology and Neoplasia 9(4), 393–404.[13] Insana, M., Liu, J., Sridhar, M. and Pellot-Barakat, C.: 2005, Ultra-sonic mechanical relaxation imaging and the material science of breastcancer, Ultrasonics Symposium, 2005 IEEE, Vol. 2, pp. 739–742.[14] Joly-Duhamel, C., Hellio, D. and Djabourov, M.: 2002, All gelatin net-works: 1. Biodiversity and physical chemistry, Langmuir 18(19), 7208–7217.[15] Kiss, M., Varghese, T. and Hall, T.: 2004, Viscoelastic characterizationof 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, PrenticeHall PTR.[18] Malvern, L. E.: 1969, Introduction to the mechanics of a continuousmedium, 1st edn, Prentice-Hall Inc., New Jersey.[19] McLaughlin, J. and Renzi, D.: 2006, Shear wave speed recovery in tran-sient elastography and supersonic imaging using propagating fronts,Inverse Problems 22(2), 681–706.105Chapter 4. Influence of Boundary Conditions[20] Nightingale, K., McAleavey, S. and Trahey, G.: 2003, Shear-wave gen-eration using acoustic radiation force: in vivo and ex vivo results, Ul-trasound in Medicine and Biology 29(12), 1715–1723.[21] Sandrin, L., Cassereau, D. and Fink, M.: 2004, The role of the couplingterm in transient elastography, J Acoust Soc Am 115(1), 73–83.[22] Sandrin, L., Tanter, M., Catheline, S. and Fink, M.: 2001, Time-resolved 2d pulsed elastography: experiments on tissue-equivalentphantoms 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 elastog-raphy, IEEE Trans Ultrasonics Ferroelectrics and Frequency Control49(4), 436–446.[24] Sarvazyan, A. P., Rudenko, O. V., Swanson, S. D., Fowlkes, J. B. andEmelianov, S. Y.: 1998, Shear wave elasticity imaging: a new ultrasonictechnology of medical diagnostics, Ultrasound in Medicine and Biology24(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 mea-sured by MR elastography, Magnetic Resonance Imaging23(2), 159–65.[26] Thompson, J. H. C.: 1933, On the theory of visco-elasticity: A thermo-dynamical treatment of visco-elasticity, and some problems of the vibra-tions of visco-elastic solids, Philosophical Transactions A 231, 339–407.[27] Turgay, E., Salcudean, S. E. and Rohling, R.: 2006, Identifying the me-chanical properties of tissue by ultrasound strain imaging, Ultrasoundin Medicine and Biology 32(2), 221–235.[28] Zahiri-Azar, R. and Salcudean, S.: 2006, Motion estimation in ultra-sound 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 andmechanical measurements of viscoelastic properties of soft tissues, Ul-trasound in Medicine and Biology 33(10), 1617–1631.106Chapter 5A High Frame RateUltrasound System for theStudy of Tissue Motions∗5.1 IntroductionMedical ultrasound imaging systems are used to image both the echogenicityof the tissue, as in B-mode imaging, and the tissue motion, as in Dopplerimaging [1,2]. Over the past two decades considerable advances have beenmade in the field of tissue motion imaging. State-of-the-art speckle trackingsoftware is capable of measuring sub-micron displacements inside tissue [3–11].A fundamental limitation of ultrasound based motion tracking is theinherent low frame rate of the ultrasound systems. Due to the finite speedof the ultrasound waves in human tissue, acquiring a full 2D B-mode imagerequirestensofmilliseconds. Alowframeratelimitsthescopeofthemotionsthat can be properly detected by the system; for a typical frame rate of100Hz, 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 cardiovascularimaging, higher frame rates are required [2]. Motions induced inside thetissue 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 theframe rate of ultrasound systems for motion estimation. In conventionalpulsed 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 ofTissue Motions.”107Chapter 5. High Frame Rate Ultrasoundframe rate, the so called pulse repetition frequency (PRF) in that context [1,2]. This frequency also determines the range of the speeds which the systemcan measure. Higher speeds demand a higher pulse repetition frequency.In conventional color (and power) Doppler imaging, the window of inter-est is subdivided into sectors [1,2]. The width of the sectors is determinedby the desired PRF, which is set by the user. Each sector is imaged repeat-edly a number of times at the PRF, and velocity estimations are carriedout. Then the system moves on to acquiring the next sector. When allthe sectors have been acquired, the velocity data for the whole window ofinterest is updated on the screen. Sonoelastography directly uses this datafor elasticity imaging [12].The idea of dividing the region of interest into a number of sectors, eachof which can be acquired separately at a higher frame rate, has been usedin a number of high frame rate imaging systems. Synchronization of theexcitation with the sector acquisition is used in [17,18] to achieve a framerate of 3kHz.Synchronization of the sector acquisition with the heart beat throughechocardiogram (ECG)-gating is used in [19,20] to achieve an effective framerate of 481Hz at an imaging depth of 11cm and 64 line density for cardiacimaging.In [21], the authors use morphing to produce additional frames betweensuccessive acquired frames. The velocity data obtained from the speckletracking is used in the morphing algorithm. Although this method helpsthe clinicians by increasing the frame rate of the video seen on the screen, itdoes not overcome the inherent Nyquist limit for measuring high frequencycomponents of the tissue motion.Another approach to high frame rate ultrasound imaging is to use non-conventional pulsing techniques. Ultrafast ultrasound is one such techniquein which transmit focusing is not applied [22]; all the probe crystals are firedat once to create a linear wave front. Dynamic receive focusing is applied toall the received signals by post processing to form the equivalent RF-linesfor each crystal. The lack of transmit focusing slightly lowers the quality ofB-mode images. However very good results have been obtained in detectinghigh 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 andprocessing of data from an aperture of crystals is required in these systems[23].A number of coded pulsing techniques have also been proposed whichcan increase the frame rate [24,25]. Frequency modulation was used to108Chapter 5. High Frame Rate Ultrasoundincrease the frame rate by a factor of 5 and Hadamard spatial encoding incombination with frequency modulation increased the frame rate by a factorof 25 [24]. A slight degradation in the B-mode image quality is reported withthese methods. However, the performance of the system in motion trackinghas not been studied.In [25] binary codes are used on the transmitted signals to produce dis-tinct impulse responses from different directions in the region of interest.The received signals are then processed in parallel to decode the data andform the image. Since multiple receive directions are compounded, blockingartifacts occur. The performance of the system in motion tracking has notbeen studied.In this article we report a high frame rate ultrasound system which usesthe concept of sector subdivision and can readily be implemented on conven-tional ultrasound systems. A detailed analysis is carried out and conditionsset forth on the tissue motions that can accurately be estimated using thissystem. The main novelty of the system is in its phase correction schemewhich makes it possible to make virtually simultaneous motion measure-ments at all points throughout the tissue being imaged. Another advantageof the system is its real-time capability. This enables the user to overlay mo-tion estimation images on B-mode images while the probe is scanning. Thesystem has been implemented on a Sonix RP machine (Ultrasonix MedicalCorporation, Richmond, BC, Canada). Examples of the applications of thesystem are given.The rest of this article is organized as follows. In Section 5.2 the highframe rate acquisition technique is reviewed and the novel phase correctionscheme presented. Section 5.3 presents a set of conditions which shouldbe satisfied by the tissue motion. These conditions ensure that the systemestimates the tissue motion correctly. Section 5.4 presents some applicationsof the system with experimental results. We conclude the article in Section5.5.5.2 Sector-Based High Frame Rate AcquisitionTechnique5.2.1 SequencingIn this section we discuss the sector-based sequencing scheme on which thehigh frame rate data acquisition is based. In this discussion, instead of theterm “firing an element” we use the term “firing a scan-line” by which we109Chapter 5. High Frame Rate Ultrasoundmean that the usual beam forming technique for focusing the beam at acertain depth is performed and an aperture of elements are fired.In conventional B-mode imaging, the scan-lines are fired one by onein 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-modeacquisition time-line for two frames on a twelve element probe. Each slantedline represents the points at different depths on one of the scan-lines. Thenumber 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 twelveelement probe.Consider a typical point at depth d on the typical scan-line No. 9. Thispoint is imaged at time t1 in the first frame and at time t2 in the secondframe. If we denote the time required to acquire a single line by ¯T, the timeinterval between successive acquisitions of the point is,t2 −t1 = Ts = 12¯T. (5.1)In motion estimation algorithms, the successive acquisitions of the RF-data from each point are used to estimate the motion at that point. Hencethe displacements are sampled at the frequency,fs = 1Ts= 112 1¯T (5.2)This frequency is the same for all the points on all the twelve scan-lines of110Chapter 5. High Frame Rate Ultrasoundthe probe in this sequencing scheme, and is equal to the frame rate of theB-mode. For a typical probe of 128 elements and a typical imaging depthof 6cm this frequency is equal to 100Hz.Figure 5.2 shows a typical sector-based sequencing with the twelve el-ement probe which would increase the Nyquist limit by a factor of 3. Inthis 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 fivein the fifth firing, the sequencer goes back to the first scan-line and repeatsthe first four scan-lines a couple of times (just once in this case for illustra-tion) before moving on to the fifth scan-line. The twelve scan-lines of theprobe are thus divided into three sectors (each, four scan-lines wide), andeach sector is scanned twice (a multitude of times in general) before the nextsector is acquired.Figure 5.2: Acquisition time-line for two frames of high frame rate acquisi-tion data on a twelve element probeAs can be seen from Fig. 5.2, the time interval between successive scan-nings of the typical point at depth d on scan-line No. 9 has been reduced. Inother words the sampling frequency has been increased during the intervalthe point is being observed (two frames in this example for illustration).111Chapter 5. High Frame Rate UltrasoundThe sampling frequency in this case is,fs = 1Ts= 14 1¯T (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 aregrouped together in a sector and the sector is scanned K times in a sequence,the sampling frequency would be,fs(M) = 1M 1¯T (5.4)From this formula, the highest frame rate is achieved with M = 1. In thiscase each line is scanned K times, before scanning the next line.Note that the total time needed to acquire K complete frames does notdepend on M, the number of scan-lines grouped together in a sector, but onN, the total number of crystals in the probe,Ttotal = parenleftbigK(M ¯T)parenrightbig NM = KN ¯T (5.5)Note that for the example of Fig. 5.2, the number of lines in each sectorM is equal to 4, the number of repetitions (frames) K is equal to 2, andthe total number of elements in the probe N is equal to 12. The number ofsectors needed to cover the entire probe is,NM =124 = 3. (5.6)As noted in the introduction, the presented sector-based sequencingscheme is used in conventional color and power Doppler imaging. In theseimaging modalities a window (region of interest) is used which reduces thetotal 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 desiredspeed range. This frequency is the same as fs(M) presented in (5.4).The presented sequencing scheme does not change the overall acquisitionspeed in terms of the total time needed for acquiring a number of frames.However it increases the sampling frequency by changing the order in whichthe lines are acquired. This reordering also changes the delay pattern be-tween the acquisition of the scan lines. In the next subsection we presentour novel phase correction algorithm for compensating these delays.112Chapter 5. High Frame Rate Ultrasound5.2.2 Time Delay CompensationThe theory of discrete signals and systems deals with discretizations ofcontinuous-time signals and systems obtained through sampling. By far,the most common type of sampling is the one in which all signals and sys-tems are sampled simultaneously and at a fixed rate. Although cases wherethe sampling rate changes or the signals are not sampled simultaneously areless common, the theory still provides the necessary tools for dealing withthem. Ultrasound based motion tracking is one.Tissue which is undergoing dynamic deformations can be considered tobe a continuous-time dynamic system. The signals of interest are the posi-tions or the velocities of the different points. It is of interest to sample allthese signals simultaneously. An ideal time-line for such an acquisition isshown in Fig. 5.3. A comparison with Fig. 5.1 shows that,• the scan-lines are no longer slanted, which means that all the pointsat 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 thescan-lines are acquired simultaneously.In reality, because of the finite propagation speed of the ultrasound insidetissue, this is impossible to achieve. The time delay from the tick of thesampling clock to the acquisition time of a point depends on the sequencingmethod used and the position of the point in the sequence.Figure 5.3: A hypothetical acquisition time-line for delay-free simultaneoussampling on a twelve element probe113Chapter 5. High Frame Rate UltrasoundFor periodic motions, the time delays can be compensated in the fre-quency domain. Assume that,• the probe does not move with respect to the tissue,• the speed of ultrasound is constant throughout the tissue and equalto c,• and that all the points to be scanned are vibrating at a single frequencyof fe. For motions containing a multitude of frequencies, the processto be presented can be repeated for each frequency separately.The acquisition process is triggered by the sampling clock at time zero. Wedenote the measured displacement phasor (at frequency fe) of the points online 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 ultrasoundhas a finite traveling speed c, the points on a single line are not all scannedat the same instant. A point at depth d inside the tissue would be scannedt = 2d/c later than the first point on the line. Since the point is vibratingat frequency fe, by this time, its phase has advanced by,φ = 2pifet = 2pife2dc . (5.8)The phasor therefore needs to be compensated as,Ucomp1(d;m,n) = Umeas(d;m,n)expparenleftbigg−j2pife2dcparenrightbigg(5.9)The next time delay to deal with is the time delay from the acquisitionof the first scan-line to the scan-line m in a sector. From Fig. 5.2, thedelay between the two lines is (m−1)¯T. Therefore, the phasor needs to becompensated as,Ucomp2(d;m,n) = Ucomp1(d;m,n)expparenleftbig−j2pife(m−1)¯Tparenrightbig (5.10)The last time delay to deal with is the time delay from the acquisitionof the first scan-line in sector one to the first scan-line in sector n. For Mlines in a sector, each scanned K times, this time delay is (n − 1)MK ¯T.114Chapter 5. High Frame Rate UltrasoundTherefore the phasor needs to be compensated as,Ucomp3(d;m,n) = Ucomp2(d;m,n)expparenleftbig−j2pife(n−1)MK ¯Tparenrightbig (5.11)The phasor Ucomp3(d;m,n) contains the amplitude and the correctedphase information for all the points, which by assumption are all vibratingat frequency fe. This phasor could be used in the frequency domain forinverse problem calculations, or converted into the time domain,u(t;m,n) = |Ucomp3(d;m,n)|cos(2pifet+∠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 16elements wide, M = 16, and 50 frames are acquired, K = 50. Assume thatthe 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,¯T = 2×77mm1540m/s = 100µ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, 2pi ×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)115Chapter 5. High Frame Rate Ultrasound5.2.3 Implementation IssuesThe scenario presented in the previous subsection assumes that the ultra-sound 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 pro-grammed into a sequencer. The sequencer then loops through that sequencerepeatedly.For the Ultrasonix RP system which is used for this research, a typicallimit of 500 scan-lines could be considered (the exact number depends ondifferent settings such as depth and sampling frequency). For a 128 elementprobe, if K = 100 samples (frames) are required, the total number of linesto 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 sequencerloop through that sector the desired number of frames, K, and then repro-gram the sequencer for the next sector. For instance if the sector size isM = 8, and K = 100 samples are required, this method requires only 8lines to be programmed and the sequencer repeats the sequence 100 times.The sequencer is then stopped and the next sector programmed. A flowchartof the process is shown in Fig. 5.4.The drawback of this solution is that there will be a (variable) timedelay between the acquisitions of the different sectors. This is due to thetime required for the software to download the new sector to the sequencerhardware. A typical time-line for this acquisition scheme is presented in Fig.5.5. The variable delay invalidates the third compensation step carried outin the previous subsection. The time delay for inter-sector compensation isnot (n−1)MK ¯T anymore, but an unknown value. Different solutions areavailable to this problem.One approach is to use time stamps. If the time that the data collectionis started for each sector is saved, say in Ti, the phasors can be compensatedas,Ucomp3(d;m,n) = Ucomp2(d;m,n)exp(−j2pife(Tn −T1)) (5.18)Another approach is to choose the sectors so that they overlap. Forinstance if the first sector consists of scan-lines 1 to 4, the second sector canconsist of scan-lines 4 to 7. By selecting the sectors in this way, the scan-line No. 4 is common to both sectors. Now there are two measurementsavailable for this line, which after the first two steps of compensation become116Chapter 5. High Frame Rate UltrasoundFigure 5.4: Flow chart of sequence programming for high frame rate imageacquisition117Chapter 5. High Frame Rate UltrasoundFigure 5.5: Acquisition time-line involving reprogramming delay on a twelveelement probeUcomp2(d;4,1) and Ucomp2(d;1,2). Since we have assumed that the probe isnot moving with respect to the tissue during the acquisition process, thesetwo phasor profiles should be exactly the same for all depths d, except for aconstant phase delay. The phase delay can be found by a simple linear leastsquares:minα,β∈CbardblUcomp2(d;1,2)−(αUcomp2(d;4,1)+β)bardbl2 . (5.19)The third compensation step becomesUcomp3(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 sectorconsists of scan-lines 1 to 4 and the second sector consists of scan-lines 5 to8. Based on scan-lines 1 to 4 a prediction can be made about the phasor ofthe points on scan-line No. 5. For instance a linear extrapolation in timefrom scan lines 3 and 4 yields,Upredicted(d;1,2) = 2Ucomp2(d;4,1)−Ucomp2(d;3,1). (5.21)118Chapter 5. High Frame Rate UltrasoundThis predicted phasor profile and the measured phasor profile after compen-sation, Ucomp2(d;1,2), should be almost the same for all depths d, exceptfor a constant phase delay. The phase delay can be found and compensatedfor by a similar procedure to (5.19) and (5.20). This is the solution we haveimplemented in our system.5.3 Requirements on the Tissue MotionAs was shown in the previous section, if all the points are vibrating at afrequency of fe in the steady state, it is possible to compensate for the ac-quisition phase lags in the phasors. Now we address the question of whetherit is possible to obtain the phasors of displacements for every frequency fe.To further elaborate on the question, consider a single point inside thetissue which is vibrating at a frequency of fe. The displacement of thispoint is sampled at a frequency of fs(M) given by (5.4). Following theNyquist criterion, to be able to recover the correct amplitude and phase ofthe displacement, it is necessary that,fe < 12fs(M). (5.22)If it were possible to observe the displacement for an infinite amount oftime and obtain infinitely many samples, (5.22) would be the only restrictionon the tissue motion. However a review of Subsection 5.2.1 reveals that foreach point, only K samples are taken, and the observation time span foreach point is,Tobs = Kfs(M). (5.23)This has very important consequences on the applicable values of fe.With reference to Fig. 5.5 each sector is being observed during a dif-ferent interval in time. Conceptually speaking, if the collected data fromthe different sectors are to be brought together and reconciled, the motionsshould be periodic with the observation interval being an integer multipleof the full period.In Fourier theory, a signal with a limited time duration Tobs is analyzedwith a Fourier series. The implicit assumptions in the conversion of such asignal 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 periodicsignal.119Chapter 5. High Frame Rate UltrasoundFigure 5.6: The extension of a limited duration signal to a periodic signalFigure 5.7: Signal processing pipeline of the high frame rate ultrasoundsystem120Chapter 5. High Frame Rate UltrasoundIn the frequency domain, the power spectrum of a periodic signal y(t)of period Tobs contains energy only at integer multiples of the fundamentalfrequency,y(t) =summationdisplayl∈Zal exp(j2pilffund), (5.24)ffund = 1Tobs. (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 limitsthe generality of the shape of the signal over the time duration Tobs. Howeverthe Nyquist criterion (5.22) limits the range of the permissible frequenciesfor the tissue motion. Substitution of (5.26) into the criterion gives−K2 < l < K2 . (5.27)In summary, the tissue motion needs to be a periodic band-limited mo-tion. Most of the natural body motions such as the motions resulting fromthe heartbeat are periodic. However the period might change slightly be-tween different pulses. A method for compensating these variations can befound in [20].If the imaging system is used together with an excitation technique forelastography applications, the total number of frequencies which could beused 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 alwaysa multiple of the fundamental frequency, the inverse of the observation pe-riod. 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 permis-sible frequencies could be used. The best combination may depend on theparticular application.For our experiments, custom software was developed that enables theuser to choose the line duration ¯T, the sector width M, and the number ofsamples K. It computes the permissible excitation frequencies and presentsthem as a list from which the user can choose the frequencies for which thephasors are computed. In between the high frame rate acquisition sequences,B-mode scans are carried out. The software can then overlay the amplitude121Chapter 5. High Frame Rate Ultrasoundor 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 atypical probe of N = 128 crystals and a typical imaging depth of d = 6cmif K = 100 samples are taken the latency is 1 second. The signal processingpipeline is shown in Fig. 5.7.5.4 Applications and Experimental ResultsIn this section applications of the high frame rate system in dynamic elas-tography are presented. In ultrasound elastography, the tissue is subjectedto a constant or time-varying force and the response of the tissue to thisforce is monitored by an ultrasound system. The response consists of thetissue motions at different locations as a function of time. The developedhigh frame rate system makes it possible to monitor tissue motions for ex-citations covering a wider range of frequencies. Moreover, since most of themotion-tracking algorithms use a similarity measure between two frames toestimate the motion, at higher frame rates (shorter intervals) the similaritiesare 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 throughwhich the ultrasound probe is fitted from below. The test sample is placedon top of the table over the ultrasound probe. The exciter consists of aspeaker and a sound tube. The sound tube acts as a wave guide betweenthe speaker and is made of plexi-glass. The inner and outer diameters are57mm and 64mm, respectively. The tube is clamped vertically to a columnand its height can be adjusted so that the open end of it touches the topsurface of the test sample. The pressure variations caused by the soundwaves inside the tube create a dynamic force excitation on the surface ofthe sample. Note the choice of the axes in Fig. 5.8; the axial, lateral andelevational directions in the ultrasound image are along the x, y and z axesrespectively. In this coordinate system, the depth d is the axial coordinatex, and the sector index variables m and n determine the lateral coordinatey:U(d;m,n) → U(x,y) (5.28)The first test is carried out on a cylindrical phantom. The phantomwas molded out of Super Soft PlasticTM (M-F Manufacturing Co., Inc. FortWorth, Tx, USA), a PolyVinyl Chloride (PVC)-based plastic. SigmaCellCellulose, Type 50 (Sigma-Aldrich Inc., St Louis, MO, USA) was added as122Chapter 5. High Frame Rate UltrasoundFigure 5.8: Schematics of (a) the side view and (b) the front view of theexperimental setup developed for testing the system. Note the choice of thecoordinate axes.123Chapter 5. High Frame Rate UltrasoundTable 5.1: Parameters of the high frame rate systemparameter valueprobe type L9–4 linear arraycenter frequency 6.5MHzsampling frequency 20MHzdepth of focus 30mmimaging depth 56mmtotal number of crystals (N) 128sector width (M) 8line duration (¯T) 100µsnumber of samples or frames (K) 50the ultrasound scatterer to develop the speckle pattern required for motiontracking. The diameter of the cylindrical phantom is 90mm and its heightis 65mm. The parameters of the high frame rate system are summarizedin Table 5.1. In this case the sampling frequency is 1.25kHz. Based on theequations (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 isused. This excitation is created using MATLAB and applied to the speakerthrough the output of the sound card on the PC. The speaker is an off-the-shelf PC-speaker and comes with an amplifier. The axial displacementsD(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 thedisplacements at the frequency of excitation is then computed,U(x,y) =k=50summationdisplayk=1expparenleftbigg−j2pik·100Hz1250HzparenrightbiggD(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 8scan-lines wide, can clearly be distinguished. Fig. 5.9(b) shows the effectof the phase compensation for the delays on each scan line, i.e. the phaseof Ucomp1(x,y). Since the compensated phase is small, the difference is notnoticeable. Fig. 5.9(c) shows the effect of the phase compensation for thedelays 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 can124Chapter 5. High Frame Rate Ultrasoundbe seen clearly to have been reduced. Fig. 5.9(d) shows the effect of thephase compensation for the delays between the sectors, i.e. the phase ofUcomp3(x,y). The phase is clearly smooth in this figure. This shows theeffectiveness 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 scan-line (c) after compensation for the delays between the scan-lines in eachsector (d) after compensation for the delays between the sectors.In the second experiment, a multiple frequency excitation is used on thesame phantom. The excitation used is of the following form:excitation = 1sin(2pi75Hzt)+1sin(2pi100Hzt)+2sin(2pi150Hzt)+3sin(2pi200Hzt) (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 canbe seen all images show smooth variations. The nodes and peaks becomemore frequent with the increase in frequency. Note that without the phasecorrection, the amplitude images would still be smooth. However, the phase125Chapter 5. High Frame Rate Ultrasoundimages would be striped as in Fig. 5.9(a).Figure 5.10: (a), (b), (c), (d) the amplitude, and (e), (f), (g), (h) correspond-ing 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 calcu-lated. The results are presented in Fig. 5.11. Note that without the phasecorrection it would not be possible to compute the lateral derivatives, ∂U/∂yand ∂2U/∂y2. Because of the delay in each scan-line, the axial derivatives∂U/∂x and ∂2U/∂x2, would also contain slight errors.The spatial derivatives of the displacement field are used extensively inelastography 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∂y + ∂V∂x , (5.31)where V(x,y) is the lateral displacement. The Laplacian of U(x,y,z) isusually approximated by:∇2U(x,y,z) = ∂2U∂x2 +∂2U∂y2 +∂2U∂z2≈ ∂2U∂x2 +∂2U∂y2 (5.32)≈ ∂2U∂x2 . (5.33)126Chapter 5. High Frame Rate UltrasoundFigure 5.11: Spatial derivatives of the axial displacement phasor U(x,y) atthe excitation frequency of 100Hz. (a) |∂U/∂x| (b) |∂2U/∂x2| (c) |∂U/∂y|(d) |∂2U/∂y2| (e) ∠∂U/∂x (f) ∠∂2U/∂x2 (g) ∠∂U/∂y (h) ∠∂2U/∂y2To demonstrate the effectiveness of the system for dynamic elastography,two inclusion phantoms were tested and their elastograms obtained. Wechose direct inversion of the wave equation for obtaining the elastograms,E(x,y) = ρ (j2pife)2U(x,y)parenleftBig∂2∂x2 +∂2∂y2parenrightBigU(x,y), (5.34)where fe, as before, is the frequency at which the phasor is computed (theexcitation frequency). Note that the factor j2pife in the numerator is thefrequency 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 itis harder than the surrounding. The blocks and the inclusions were madefrom Super Soft Plastic with different concentrations of Plastic Softener orPlastic Hardener. A schematic of the phantom geometry and imaging planeis presented in Fig. 5.12.The setup of Fig. 5.8 is used to excite and image the phantoms. The127Chapter 5. High Frame Rate UltrasoundFigure 5.12: Outline of the geometry of the inclusion phantoms used in theelastacity experiments.excitation used is of the following form:excitation = 1sin(2pi100Hzt)+1.5sin(2pi125Hzt)−1.5sin(2pi150Hzt)−2sin(2pi175Hzt)+2sin(2pi200Hzt). (5.35)For each of the excitation frequencies, the phasor U(x,y) was calculatedand was used in (5.34) to obtain an elastogram E(x,y). The obtained elas-tograms for different frequencies were then averaged to obtain the final elas-togram.The results for the hard and soft inclusion phantoms are presented inFig. 5.13 and 5.14 respectively. Since the echogenicity of the inclusions areclose to their surrounding material, it is difficult to recognize them in the B-mode images. However, they can clearly be recognized from the elastogramsin both cases.5.5 Summary and ConclusionA high frame rate ultrasound system is introduced in this article whichcan readily be implemented on conventional ultrasound systems without128Chapter 5. High Frame Rate UltrasoundFigure 5.13: (a) B-mode and (b) elastogram images of a phantom with acircular hard inclusion.Figure 5.14: (a) B-mode and (b) elastogram images of a phantom with acircular soft inclusion.129Chapter 5. High Frame Rate Ultrasoundthe need for any hardware modifications. The system uses the concept ofsector 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 tobe scanned.A phase correction scheme was presented to compensate for the timedelays involved in the acquisition process. Necessary conditions on the fre-quency content of the tissue motion were presented, in order for the dis-placements to be recovered properly.Since the system simply reorders the pulse sequence, the tissue is stillimaged one line at a time. As a result the total acquisition time with thissystem can be a fraction of a second. The probe and tissue should be sta-tionary during this acquisition interval to prevent motion artifacts. Thesystem can be used together with periodic excitations to study tissue char-acteristics, and as such is not meant as a replacement for other high speedimaging techniques used in applications with different requirements.The phase compensation scheme presented here can readily be extendedto the 3D case. For 2D phased array transducers the extension is straightforward. For mechanically actuated 3D probes with linear phased arraytransducers, the motion of the stepper motor inside the probe can be syn-chronized with the sector acquisition to produce known delays between dif-ferent image planes. These known delays can then be compensated for in asimilar manner.The high frame rate system can be used together with the exciters usedin elastography techniques, such as mechanical exciters and acoustic radi-ation force (ARF) exciters, to obtain displacement data of the tissue atdifferent frequencies of excitation. The excitation frequency can go up toa few kiloHertz and the system can still estimate the displacement data.These data can then be used in the inverse problem solving techniques torecover the tissue elasticity and viscosity.130References[1] R. Cobbold, Foundations of Biomedical Ultrasound. NY: Oxford Uni-versity Press, 2007.[2] J. Jensen, Estimation of Blood Velocities Using Ultrasound: A SignalProcessing Approach. Cambridge University Press, 1996.[3] I. Hein and W. O’Brien, “Current time-domain methods for assessingtissue 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 ofnormal and shear components of the 3d strain tensor in elastography,”Phys. Med. Biol., vol. 45, pp. 1553–1563, 2000.[5] R.Righetti, J.Ophir, andP.Ktonas, “Axialresolutioninelastography,”Ultrasound Med. Biol., vol. 28, no. 1, pp. 101–113, 2002.[6] S. Langeland, J. Dapos, H. Torp, B. Bijnens, and P. Suetens, “A simu-lation study on the performance of different estimators for two dimen-sional velocity estimation,” in Proc. IEEE Ultrasonics Symp., vol. 2,2002, p. 18591862.[7] R. Righetti, S. Srinivasan, and J. Ophir, “Lateral resolution in elastog-raphy,” Ultrasound Med. Biol., vol. 29, no. 5, pp. 695–7–4, 2003.[8] F. Viola and W. Walker, “A comparison of the performance of timedelay estimators in medical ultrasound,” IEEE Trans. Ultrason., Fer-roelect., Freq. Contr., vol. 50, p. 392401, 2003.[9] R. Zahiri and S. Salcudean, “Motion estimation in ultrasound imagesusing time-domain cross-correlation with prior estimates,” IEEE Trans.Biomed. Eng., vol. 53, no. 10, pp. 1990–2000, 2006.131Chapter 5. High Frame Rate Ultrasound[10] G. Pinton, J. Dahl, and G. Trahey, “Rapid tracking of small dis-placements 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 ultra-sound echo signals using individual sample tracking,” IEEE Trans. Ul-trason., Ferroelect., Freq. Contr., vol. 55, no. 12, pp. 2640–265, 2008.[12] L. Taylor, B. Porter, D. Rubens, and K. Parker, “Three-dimensionalsonoelastography: Principles and practices,” Phys. Med. Biol., vol. 45,pp. 1477–1494, 2000.[13] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulusimaging with 2-d transient elastography,” IEEE Trans. Ultrason., Fer-roelect., Freq. Contr., vol. 49, no. 4, pp. 426–435, April 2002.[14] J. Bercoff, R. Sinkus, M. Tanter, and M. Fink, “Supersonic shear imag-ing: 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 inimages generated with transient acoustic radiation force,” Ultrasoundin Med. & Biol., vol. 32, no. 1, p. 6172, 2006.[16] J. J. Dahl, G. F. Pinton, M. L. Palmeri, V. Agrawal, K. R. Nightin-gale, and G. E. Trahey, “A parallel tracking method for acoustic radia-tion 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 time-varying mechanical viscoelastic parameters of mimicking deep veinthrombi with 2d dynamic elastography,” in Proc. of IEEE Ultrason-ics Symposium, 2007, pp. 1009–1012.[18] A. H. Henni, C. Schmitt, and G. Cloutier, “Analytical modeling ofplane shear wave diffraction by a radially layered cylinder for dynamicvascular 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 imag-ing technique for high frame-rate and full-view cardiovascular ultra-sound and elasticity imaging,” in Proc. of IEEE Ultrasonics Sympo-sium, 2007, pp. 880–883.132Chapter 5. High Frame Rate Ultrasound[20] S. Wang, W.-N. Lee, J. Provost, J. Luo, and E. E. Konofagou, “Acomposite 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 framerate 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 IEEEUltrasonics Symposium, 2002, pp. 1811–1820.[23] C. Fabian, K. Ballu, J. Hossack, T. Blalock, and W. Walker, “Devel-opment 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 inmedical 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 forhigh-speed ultrasonic imaging,” IEEE Trans. Med. Imaging, vol. 17,no. 6, pp. 923–934, 1998.133Chapter 6Conclusion6.1 SummaryThe inverse problem of elastography was studied in this research [1–14].The objective of the inverse problem is to find the elasticity of the tissueas a function of position, given the displacements as a function of positionand time. In ultrasound elastography, usually only one component of thedisplacement, namely the axial component, is available from the measure-ments. Moreover, the measurement is restricted to the plane of imaging. Inthis case, the inverse problem is that of finding tissue elasticity based onpartial displacement measurements.The wave equation approach to the inverse problem, in the case of partialdisplacement measurements, addresses the problem by finding the phasespeed of the waves [2–5,7,8,10,11,13,14]. In this thesis, the limitationsof this approach were clearly defined. In particular, it was shown thatcontrary to common assumptions, longitudinally polarized waves may notalways travel with the so called shear wave speed in a medium. The shearwave speed is the phase speed of plane shear waves in an infinite mediumand is given bycs =radicalbiggµρ . (6.1)It was shown that in a finite block of tissue, when a large profile mechanicalexciter is used to compress the tissue longitudinally, the phase speed is givenbycphase =radicalbigg3µρ , (6.2)which is the speed of longitudinal waves in thin rods.A further analysis of the phase speed revealed that the phase speed, inaddition to the tissue elasticity, depends on the shape of the wave beam, andon the boundary conditions of the medium. Therefore a universal relation-ship between the phase speed and the tissue elasticity cannot be assumed.However in particular cases and under certain boundary conditions, the re-134Chapter 6. Conclusionlationship (6.1) holds and the wave speed can be used to estimate tissueelasticity.Different methods have been used for estimating the phase speed. Thesemethods are based on local frequency estimation. For instance they measurethe wavelength of the wave inside the medium. For accurate measurementof the wavelength, it is of interest to have a few periods of the wave observedin the length of the region of interest. To place more periods in the region ofinterest, the frequency of excitation should be increased. This is a commonfeature of all phase estimation methods; to achieve a higher resolution, ahigher frequency of excitation should be used.In this thesis a technique was presented for increasing the frame rateof an ultrasound system to image higher frequencies of excitation. Theconventional method of ultrasound imaging has a limited frame rate whichdepends on the depth of imaging, line density, and width of the regionof interest. Typical frame rates are in the range of 50Hz to 200Hz. Bythe Nyquist criterion, this limits the frequency of excitation to 25Hz to100Hz. The presented method in this thesis was easily implemented on aconventional ultrasound system, and increased the frame rate to a range of1kHz to 10kHz.The presented high frame rate system was used to successfully trackmotions with a frequency of up to 500Hz. The tracked displacements wereused to find the phase speed, and elasticity. The elastograms obtained fromhard inclusion and soft inclusion phantoms clearly revealed the inclusion.At this point, it can be stated that all the primary objectives of thisthesis were achieved.• A method was devised for accurately measuring the elasticity of tissueand tissue mimicking material samples for setting gold standards.• Tissue mimicking materials were fabricated with controllable elastic-ity.• The relationship between the phase speed and tissue elasticity wasestablished through a theoretical analysis of the elastic wave equation.• An ultrasound elastography system was developed which uses the de-veloped inversion algorithms for forming elastograms.Therefore the hypothesis that a wave equation based approached can be usedin ultrasound elastography for forming elastograms has been confirmed.135Chapter 6. Conclusion6.2 Contributions6.2.1 RheometryA novel framework was devised in Chapter 2 for viscoelastic characteriza-tion of tissue and tissue mimicking material. Longitudinal wave excitationwas used on blocks of test phantoms and displacements were imaged byan ultrasound system. Two methods were presented for analyzing the dis-placements. The first method estimated the Young’s modulus based on thedistance from the peak to the node of the standing wave patterns inside theblock. The second method used a model fitting scheme to find the values ofthe Young’s modulus and the relaxation time. The results were comparedto conventional rheometry results obtained using force-displacement mea-surements. The error in the estimated elasticities is lowest for the range of30-80kPa. The maximum error for the peak to node method was 11% and forthe model fitting method was 20%. For the viscosities, only the trends couldbe compared as the measurements were over different frequency ranges.6.2.2 Poly Vinyl Chloride PhantomsIn Chapter 2 fabrication of homogeneous blocks of tissue mimicking mate-rial was reported from Super Soft PlasticTM, a Poly Vinyl Chloride basedphantom using different proportions of the softener and hardener. SuperSoft PlasticTM was shown to be a promising phantom material for elastog-raphy studies, as its Young’s (and hence shear) modulus can be adjustedover a wide range (15–150kPa), without changing its density and the speedof ultrasound cinf. Figure 6.1 summarizes the Young’s modulus obtained vs.the volume portion of hardener or softener used in the fabrication processfor this material. However, the relaxation time and its frequency dependentbehavior is not dependent on the volume portion of Softener or Hardener(refer to Fig. 2.5(b) in Chapter 2). Therefore without any additives, thismaterial by itself is unsuitable for creating phantoms with significant vis-cosity contrast.6.2.3 Phase Speed and Tissue ElasticityA theoretical analysis of the phase speed and its relationship to tissue elas-ticity was carried out in Chapter 3 which enhanced our understanding ofthe wave phenomena. From this analysis, it was concluded that the prop-agation of elastic waves trapped in finite media, such as soft tissue, is acomplex phenomenon even without considering the viscous and nonlinear136Chapter 6. ConclusionFigure 6.1: Obtained Young’s modulus as a function of the volume portion ofPlastic Softener or Hardener added during the fabrication of the phantomseffects. The term “wave” is generally applied to any vibration pattern in-side the tissue. However the phase velocity cannot be defined for all thecases of the vibration patterns. Usually multiple waves traveling in differ-ent directions superimpose to create the vibration pattern. In this case thephase velocity might be well-defined for each wave, but not for the resultantof the superposition.Even in the simple case where a single frequency wave is traveling in asingle direction, the geometry of the wave and the medium create multiplepermissible phase speeds, and thus make it difficult to recover a meaningfulphase speed for the traveling wave. Therefore no universal relationship canbe proposed which would relate the phase speed to tissue elasticity.Based on the analysis carried out the following measures were recom-mended to minimize the artifacts in elastograms obtained using phase speedestimation. The main goal should be to create a single wave propagating inone direction at a single speed:• The excitation amplitude should be chosen small enough to reduce thereflected wave amplitudes to the level of other measurement noises.The damping present in the tissue helps by reducing the amplitude ofthe reflected waves from the boundaries of the tissue and bones.• The frequency of excitation should be kept to a minimum to preventthe appearance of higher modes. However lower frequencies result inlower damping, and therefore a trade off exists here.• The excitation pattern should be chosen so as to excite mainly the137Chapter 6. Conclusionlowest mode.6.2.4 High Frame Rate Ultrasound SystemA novel high frame rate ultrasound system was introduced in Chapter 5which could readily be implemented on conventional ultrasound systemswithout the need for any hardware modifications. The system uses theconcept of sector subdivision for sequencing which is widely used in Dopplerimaging to achieve high frame rates. For a sector of width one, i.e. a singlescan-line, the frame rate is only limited by the time it takes for the deepestpoint 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 compen-sate for the time delays involved in the acquisition process. Necessary condi-tions on the frequency content of the tissue motion were presented in orderfor the displacements to be recovered properly.Since the devised system simply reorders the pulse sequence, the tissueis still imaged one line at a time. As a result, the total acquisition timewith this system can be a fraction of a second. The probe and tissue shouldbe stationary during this acquisition interval to prevent motion artifacts.The system can be used together with periodic excitations to study tissuecharacteristics, and as such is not meant as a replacement for other highspeed imaging techniques which are used in applications with different re-quirements.6.2.5 Ultrasound Elastography SystemThe high frame rate ultrasound system was used to devise an ultrasoundelastography system, which was reported in Chapter 5. The system usesperiodic excitations, with multiple excitations frequencies to vibrate thetissue. Thehighframerateultrasoundsystemimagesthetissueandprovidesthephasecompensatedphasordisplacementsatthefrequenciesofexcitation.The phasor displacements are then used in an inversion algorithm, namelythe algebraic inversion of the wave equation, to estimate the phase speed andhence tissue elasticity. The results obtained from different frequencies areaveraged to give the final elastograms. The performance of the system wastested on hard inclusion and soft inclusion phantoms. The images obtainedclearly show the geometry and location of the inclusions through changesin relative elasticity. However, measuring absolute values of elasticity hasbeen a challenge which has not been fully addressed in this thesis. Theunderlying problem is that of estimating the local frequency (or wavelength)138Chapter 6. Conclusionof a function from a few samples of it at each of these samples. Using higherfrequencies of excitation helps, but it also presents two new challenges. Asthe frequency goes higher, the damping of the wave becomes more severe inreal tissue. Frequencies of above 500Hz cannot penetrate very deep insidereal tissue. The second problem is the appearance of multiple modes andphase speeds. This problem could be addressed by using the curl of thedisplacements, as suggested in the next section.6.3 Future WorkIn this thesis, a high-frame-rate elastography system was reported. The sys-tem works with periodic excitations in the range of 1-500Hz. The elasticityimage update rate depends on the width of the ROI and the number offrames acquired and can be as low as a few images per second. The perfor-mance of the system was tested on tissue mimicking material. Along thisline of research, a number of directions can be pursued:• Different excitation and inversion methods can be tested together withthe high-frame-rate system to find an optimal combination for per-forming in vivo tests.• The parameters involved, such as the number and values of the fre-quencies of excitation, the frame rate of the high-frame-rate system,and the number of frames to be acquired could be optimized based onthe 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 thesignals is an inherently difficult problem. This is the problem of phasespeed estimation for a small region of interest composed of for instance 100pixels by 100 pixels. It is therefore of interest to circumvent the phase speedestimation, and use the wave equation directly to find the values of themechanical properties. Along this line of research, a number of directionscan be pursued:• The wave equation can be written for the shear component of the waveby application of the curl operator on the displacement data. To beable to use the curl operator, displacements should be measured inat least two directions and over a volume. The current research atthe Robotics and Control Lab is addressing the problem of two direc-tions. By using beam steering in conjunction with the high-frame-rateand phase compensation scheme of this thesis, lateral motions can be139Chapter 6. Conclusionmeasured accurately in addition to axial motions.• Acquisition of a 3D volume of beam steered ultrasound data with axialand lateral motion tracking using 3D ultrasound probes is a furtherstep.• Different wave equation inversion algorithms based on the curl of thedisplacements can be implemented and their performance assessed ontissue mimicking material and real tissue.• Measuring absolute values of tissue elasticity is another challengingarea which requires further investigation.140References[1] R. Sinkus, J. Lorenzen, D. S. amd M. Lorenzen, M. Dargatz, andD. Holz, “High-resolution tensor mr elastography for breast tumor de-tection,” Physics Med. Biol., vol. 45, pp. 1649–1664, 2000.[2] A. Manduca, T. Oliphant, M. Dresner, J. Mahowald, S. Kruse, E. Am-romin, J. Felmlee, J. Greenleaf, and R. Ehman, “Magnetic resonanceelastography: Non-invasive mapping of tissue elasticity,” Medical ImageAnalysis, vol. 5, pp. 237–2540, 2001.[3] T. Oliphant, A. Manduca, R. Ehman, and J. Greenleaf, “Complex-valued stiffness reconstruction by for magnetic resonance elastographyby algebraic inversion of the differential equation,” Magnetic Resonancein Medicine, vol. 45, pp. 299–310, 2001.[4] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Shear modulusimaging with 2-d transient elastography,” IEEE Trans. Ultrason. Fer-roelect. 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 us-ing 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 andelastic properties ofsoft tissue using supersonic shear imaging,” in IEEEUltrasonics Symposium, 2003, pp. 925–928.[8] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity andviscosity from measurement of shear wave speed dispersion,” J. Acoust.Soc. Am., vol. 115, no. 6, pp. 2781–2785, June 2004.141Chapter 6. Conclusion[9] J. Bercoff, R. Sinkus, M. Tanter, and M. FinkJ, “3d ultrasound-baseddynamic and transient elastography: First in vitro results,” in IEEEUltrasonics Symposium, 2004, pp. 28–31.[10] J. Bercoff, R. Sinkus, M. Tanter, and M. Fink, “Supersonic shear imag-ing: 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. Abouelka-ram, and J. Culiolic, “Measurement of viscoelastic properties of ho-mogeneous soft solid using transient elastography: An inverse problemapproach,” J. Acoust. Soc. Am., vol. 116, no. 6, pp. 3734–3741, Dec2004.[12] J. Fromageau, J.-L. Gennisson, C. Schmitt, R. L. Maurice, R. Mon-grain, and G. Cloutier, “Estimation of polyvinyl alcohol cryogel me-chanical properties with four ultrasound elastography methods andcomparison with gold standard testing,” IEEE Trans. Ultrason. Fer-roelect. Freq. Contr., vol. 54, no. 3, pp. 498–509, 2007.[13] J.-L. Gennisson, S. Catheline, S. Chaffa¨ı, and M. Fink, “Transient elas-tography in anisotropic medium: Application to the measurement ofslow 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 methodsfor the detection and localization of inclusions in structures utilizingdynamic surface displacements,” in Proceedings of SPIE Vol. 5503, vol.5503, 2004, pp. 367–374.142Appendix AConventional RheometryThe measurement of the viscoelastic properties of the tissue phantoms at thefrequencies considered in this article (between 50Hz to 100Hz) by conven-tional rheometry (force-displacement measurements) is a nontrivial problem.Here we elaborate on this issue.In conventional compression based force-displacement rheometry, a blockof test material with known dimensions w × l × h is placed between twoparallel plates (see Fig. A.1). One of the plates is actuated and can moveto compress the block. The force F needed to compress the material andthe displacement of the moving plate d are measured. The stress and strainare found as,epsilon1xx = dh , (A.1)τxx = Fw×l . (A.2)The elasticity (stiffness) of the material is found as,E = τxxepsilon1xx. (A.3)Conventional rheometry can also be used to study the behavior of thestiffness as a function of frequency. For this purpose the actuated plate ofthe device is excited by a time varying waveform, causing a time varyingstrain and a time varying stress. In this case, the stiffness would beS(jω) = F{τxx(t)}F{epsilon1xx(t)}. (A.4)where F() is the Fourier transform. The extension of static rheometry to adynamic setting raises a number of very important issues which affect themeasurement results. The first of these issues is the mounting of the forcesensor. In the static case, the force sensor can be fitted between either ofthe plates and the rest of the device. In equilibrium the forces are different143Appendix A. Conventional Rheometryonly by the weight of the test block and the weight of some parts of the forcesensor. This difference can be measured and the device can be calibratedfor this setting. In the dynamic case, however, there are significant issuesrelated to each of the possible fittings.Assume that the force sensor is fitted between the actuator of the deviceand the moving plate (Fig. A.1). In this case the measured force by thesensor not only includes the applied force on the material block but alsothe inertial forces necessary to accelerate the plate and the mass of theblock. The device needs to be calibrated, but the calibration is much morecomplicated than in the static setting. The inaccuracies involved in themeasurements make this method very challenging to implement.Figure A.1: Simple schematic of the rheometry device in a conventionalrheometry configuration. In this configuration the force sensor measuresthe force on the moving plate and moves with it.Another option is to fit the force sensor between the fixed plate and thechassis of the device (Fig. A.2). Although in this case the measured forceis the force applied to the block to compress it from above, a second factorcomes into play which complicates the problem. Due to the finite speed ofpropagation of longitudinal waves in the block, there would be a time delay(and hence a phase difference) between the measured force on the fixed endof the block and the force on the moving end of the block. The higher thefrequency of the excitation, the larger the phase difference:φ(ω) = ωt = ω hcfin. (A.5)If the velocity of the longitudinal waves cfin were known, this formula could144Appendix A. Conventional Rheometrybe used directly to compensate for the phase difference,S(jω) = F{τxx(t)}F{epsilon1xx(t)}exp(jφ(ω)), (A.6)but this velocity is itself a function of the stiffness (2.30):cfin =radicalBiggEρ =radicalBiggreal(S(jω))ρ . (A.7)This makes the measurement more challenging, but still feasible for lowfrequencies.Figure A.2: Simple schematic of the rheometry device in a conventionalrheometry configuration. In this configuration the force sensor measuresthe force on the static plate and does not move.At higherfrequencies, the problem becomes more severe. Infactequation(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 volumeare perpendicular forces to the top and bottom surfaces, the strain and stressin the vertical direction are proportional, the proportionality factor beingthe Young’s Modulus E. Whether this equation could be applied to a blockof material depends on how similar the conditions are to the infinitesimalsetting.As far as the frequency is concerned, in the infinitesimal setting, theheight of the block is always negligible with respect to the wave length ofa propagating wave. This is not always true for a block. In other words,the block cannot always be considered as a lumped element with lumpedmechanical properties equal to that of its infinitesimal elements. If the wave145Appendix A. Conventional RheometryFigure A.3: The infinitesimal volume which is used in the definition of Elas-tic modulus.length is comparable to the height of the block, equation (A.4) would nolonger be valid for the block.One way to deal with the problem of the block not being a lumped ele-ment at higher frequencies is to reduce the height of the block, h. In otherwords, thin slabs could be used for rheometry measurements. But goingback to the question of the similarity between the block and the infinites-imal settings, the boundary condition that “the only forces acting on theinfinitesimal volume be perpendicular forces to the top and bottom sur-faces” can never be achieved with the block; without any lubricant on thetop and bottom plates of the device, the material would stick to the platesbecause of the friction and there would be tangential forces acting on thetop and bottom surfaces of the block. The middle portion of the block,however, resembles the infinitesimal setting to a good extent. But reducingthe block height h reduces the effective middle portion. Putting some lu-bricant helps but never removes the friction. It might be worth mentioningthat this problem mainly appears for nearly incompressible material such asliving tissue and tissue phantoms. These materials have a high tendency toexpand laterally as they are compressed.From this discussion, it follows that rheometry of soft tissue by con-ventional compression type force-displacement measurement at higher fre-quencies is not a trivial problem. In fact the proposed ultrasound basedrheometry technique can be considered as a simple solution for this prob-lem. Most commercially available rheometry devices use shear rheometrybut their bandwidth is limited to 50Hz or bellow. Even with higher band-width systems, similar issues would exist with soft tissue since the test blockswould not be lumped at higher frequencies.In this report we validated our estimated elasticity values with the values146Appendix A. Conventional Rheometryobtained from low frequency rheometry using the second configuration (Fig.A.2). The relaxation times however could not be validated directly, exceptfor the general trends.147Appendix BTissue Mimicking Materialsfor Ultrasound ElastographyIn this appendix, a brief review of some of the most commonly used tissuemimicking materials is presented. Also a detailed recipe for fabricating PVCphantoms is given.B.1 Literature ReviewA comprehensive study of hydrogels, namely gelatine and agar based phan-toms has been carried out in [1]. The purpose of the study is assessmentof the acoustic and mechanical properties of these phantoms for ultrasoundelastography. The authors mainly study Young’s modulus and sound speed.They give exact recipes for preparing the phantoms, and carry out differentmeasurements and report the results. The phantoms are based on gelatineor agar solution in a mixture of dionized water and n-propanol. n-propanolis used to increase the sound speed. Graphite or glass micro-beads are usedas scatterers. Some of the methods they use in preparation of the phantoms,such as use of vacuum chambers for degasing the mixtures, or rotating themixtures to keep the graphite afloat while the gel congeals are not reallyappropriate for making complex anatomical phantoms, as they might be forcylinderical or cubic phantoms. The authors also study the effect of timepassage on the degradation of the prepared phantoms.Another study on phantoms appropriate for ultrasonic elastography iscarried out by Madsen et.al. [2]. The authors perform an informative liter-ature survey on the state of the art phantoms used in ultrasonography andelastography. They note that most of the these studies lack a tracking of theaging of the phantoms and the changes in their characteristics. As a resultthe phantoms are not stable during the development of imaging techniquesand hence these researches suffer a methodical problem. They also proposegelatine based phantoms with dispersed saﬄower oil for ultrasonic elastog-raphy. Exact recipes are given. Different proportions of ingredients are used148Appendix B. Tissue Mimicking Materialsto achieve different characteristics. The group studies the stability of acous-tic and mechanical properties of the phantoms, over an extended period oftime in [3]. They also study the stability issue for agar and gelatine basedphantoms without oil dispersions [4].A comprehensive literature survey has been carried out by Hassan andPeppas on Poly(vinyl alcohol) polymers in [5]. This material has a widevariety of applications in biomedical engineering from contact lenses to ar-tificial implants, such as artificial cartilages. Therefore numerous methodshave been developed for specific applications. The paper does not aim togive a certain recipe for preparation of Poly(vinyl alcohol) gels, rather givesan overview of different methods and material used for obtaining differentgels.Another seminal paper has recently been published by Madsen et.al. [6].In this paper the authors present the fabrication of two anthropomorphicbreast phantoms for Ultrasound and MRI Elastography. Six tissue typeshave been prepared. The material used is gelatine with dispersed saﬄoweroil. A well-detailed recipe is given as well as the involved molding procedureused to prepare the phantom. The paper sets a good standard for preparinganthtopomorphic phantoms. The only downside of the fabricated phantomis that it should be kept in saﬄower oil, even during experiments to preventdesiccation.In [7] the authors present a bone-mimicking material base on plaster,CMC vegetable adhesive and chopped flax which shows attenuation andspeed characteristics similar to bovine bone.Another bone mimicking tissue has been proposed in [8] using apatiteand PVA. Sound speed, characteristic impedance and the density of thephantom matches those of the bovine cancellous bone.In [9] the authors use a 13% gelatine, 10% graphite, 10% formaldehydephantom with an embedded cellophane locator marker and vibrate it tocreate strain-rate (SR) images using doppler-ultrasound. The phantom isbased on the work of Hall et al. in [1]. The claim is that the mixture closelysimulates the stiffness and scattering properties of soft tissue.The authors of [10] investigate the temperature and frequency depen-dence of four types of commercially available phantoms: urethane rub-ber from ATS Labs (Bridgeport, CT), two types of condensed milk fromGammex-RMI (Middleton, WI) and ZerdineTM from CIRS Inc. (Norfolk,VA). They also study an agar based phantom developed through a Euro-peanCommission(EC)-fundedproject. TheymentionthattheInternationalElectrotechnical Commission 1390 (IEC 1996) and American Institute ofUltrasound in Medicine Standard 1990 (AIUM 1990) standards for TMMs149Appendix B. Tissue Mimicking Materialsrecommend an acoustic velocity of 1540 m/s, attenuation coefficients of 0.5dB/(cm MHz) and 0.7 dB/(cm MHz) for the frequency range 2 to 15 MHzwith a linear response of attenuation to frequency, f1.PolyVinyl Alcohol Cryogel (PVA-C) phantoms which have undergonedifferent numbers of freeze and thaw cycles to adjust their acoustic andmechanical properties are used for vessel elastography in [11]. The authorscite [12] as their reference for building the phantoms, and do not mentionedthe exact method they have used, except the number of freeze and thawcycles.In [13] PVA based phantoms are proposed which have optical, acousticand mechanical properties similar to breast tissue. The authors are perform-ing Ultrasound Assisted Optical Elastography (UAOE), a new approach toelastography. They describe in detail the recipe for making the PVA phan-toms and the freeze-n-thaw cycles used. Their work is based on [5] and [14].They also describe methods for measuring optical, acoustic and mechanicalproperties of the phantoms. The measured properties for two different PVAsand 5 different numbers of freeze-n-thaw cycles are given.Simple PVA based material are used to fabricate a number of antropo-morphic phantoms in [15]. They have used usuall freeze-n-thaw cycles toincrease the stiffness of their gels. They also study the effect of adding colorto dye the phantom on some of the characteristics.Complex tissue mimicking phantoms are proposed in [16] for multimodalimaging of prostate. These contain water, agarose, lipid particles, protein,Cu++, EDTA, glass beads, and thimerosal (preservative). Nevertheless themechanical properties of the phantoms are not of concern to the authors.Nonlinear mechanical properties of an agar-based and a gelatine-basedphantom under large static strains are reported in [17]. The authors usedbleach 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 assesmentof tissue elasticity and viscosity within the framework of shear wave elasticityimaging and acoustic remote palpation. Glycerol was added to gelatine toincrease the viscosity of the phantoms. Shear wave speed, speed of soundand 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 isreported in [19]. The work is of particular interest, as it describes the meth-ods used for integration of the anatomical features as well as the recipe forpreparing them. Different tissue-mimicking parts are fabricated using differ-ent consentrasions of Al2O3 and SiC particles inside an agar-based hydrogel.150Appendix B. Tissue Mimicking MaterialsA 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 elas-tography phantoms is polyacrylamide [20]. The authors have used TiO2to adjust the acoustic properties of the phantom, and different concentra-tions of acrylamide to alter the elasticity. The recipe for using this materialis given. The negative point about this material is that although polyacry-lamide is nontoxic, but since the polymerization is not perfect, some residualacrylamide is always present, and this latter is a neurotoxin. Nevertheless itseems that small quantities of acrylamide is present in fried food and wetherit is a carcinogen is still a point of debate. Special handling practices arerecommended.Another PVA based phantom for photoacoustic mammography has beendeveloped in [14]. In this imaging technique, the ultrasonic waves are pro-duced 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 de-sired photoacoustic properties. In a second method they use a mixture ofwater and dimethylsulphoxide to dissolve the PVA, and reach at a transpar-ent gel. Then by addition of borosilicate glass microspheres they alter theturbidity of the phantom. Except for the densities no study of mechanicalproperties is carried out.Kondo and Fujimoto have reported oil-gel phantoms [21]. The mainadvantage of these phantoms is that they do not rust. The phantoms arebased on ethylene glycol (car radiator antifrost) and propylene glycol orpolypropylene 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 sexuallubricants. The author only report sound speed, density and some relativemeasure of hardness. No recipes are given.Madsen et.al. propose a phantom geometry for testing Ultrasound andNMR elastography methods in [22]. The phantom consists of a series ofspherical balls of different sizes embedded in a background. Three differentphantoms are fabricated with this geometry. Two of them are agar-gelatinebased and the third one involves dispersions of different concentrations ofmicroscopic oil droplets in a gelatin matrix. The exact percentages of thematerial are given, but for the recipe, the authors refer the reader to [2]and [4]. The paper mainly concentrates on the molding procedure used toobtain the phantoms.Polyvinyl chloride(PVC) is another polymer which has the potential ofbeing used as a base for elastography and ultrasonography phantoms. To thebest of our knowledge, this material has not been investigated thoroughly151Appendix B. Tissue Mimicking Materialsin the literature. Although some reports on the use of this material exist[23] [24] [25] and [26].B.2 PVC Phantom RecipeB.2.1 IntroductionSuper Soft PlasticTM phantoms are made by mixing appropriate propor-tions of Super Soft PlasticTM with either Plastic SoftenerTM or PlasticTM .The ultrasound speed is pretty much independent of the mixture ratio andis around 1420m/s, slightly lower than the ultrasound speed in real tissuewhich is 1540m/s. The amount of softener or hardener however determinesthe elasticity or stiffness of the final phantom. Without adding any par-ticles, the phantoms would be transparent to ultrasound. Cellulose withaverage 50micron particle size is added to create scatterer pattern. Theratio of scatterer to plastic volume should not exceed 4%, otherwise exces-sive scattering would limit the penetration depth of the ultrasound to a fewcentimeters. Good results are obtained using 1% cellulose which is about1/4tsp for 100mL plastic.B.2.2 Preparation StepsThe following steps should be taken.• The plastic, softener and hardener container should be shaken wellbefore pouring the plastic, since the plasticizer and the liquid plastictend to form a two-phase mixture.• The plastic, softener or hardener is measured and poured inside abeaker.• The beaker is heated over the hotplate while trying to maintain auniform temperature distribution in the liquid by mixing.• When the plastic turns into a transparent liquid, the cellulose is addedslowly while stirring to prevent creation of lumps.• Let the liquid stay at a high enough temperature to keep it in theliquid state to allow for the air bubbles to surface. This may take 10to 30 minutes.• Pour the liquid slowly into the mold to allow for the generated airbubbles during the pouring process to remain on the surface and notget buried.• Let the mold to cool down well before trying to remove the phantom.Remove the phantom slowly by trying to disconnect its surfaces from152Appendix B. Tissue Mimicking Materialsthe mold surfaces.B.2.3 Other Important Notes and TipsThe following notes have been made from several attempts at making PVCphantoms:• Since the liquid plastic is heated close to its evaporation point, a lot offume is generated which may pose health and safety risks to you andthe other people working in the room. The whole process of heating,molding and cooling of the phantoms should therefore be performedunder a fume hood. The company providing the material sells themfor making fish lures. They recommend fabrication under the kitchenfume hood. I would recommend using a chemical lab fume hood forcomplete protection.• There is no water-based solvent for PVC. In other words commercialsoaps, 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. Ace-tone is very evanescent however and its vapor also poses health andsafety risks. The surface to be cleaned should be wiped with acetone.• Hotplates come in a number of different types. The ones with control-lable plate temperature and an even temperature distribution tendto be the best. However since these hotplates are not available inevery case, what happens is that the plate temperature usually goesabove the necessary temperature for heating the plastic. If the beakeris directly placed over the hotplate, this may cause the liquid at thebottom of beaker to burn and stick to the beaker bottom and alsodarkens the color of the phantom. To get a milder temperature dis-tribution in this case a pan can be placed between the plate and thebeaker. A better method would be filling a pot with sand and dippingthe beaker inside the sand and heating the pot. This would create amore even temperature distribution inside the beaker.• When retaining the plastic in the molten state to let the air bubbles tosurface, you don’t want the temperature to go very high, because thelighter chemicals would evaporate and change the properties of yourphantom. Also the phantom would start to darken too much and looseits visual appeal. On the other hand too low a temperature wouldcreate a very viscose environment in which small air bubbles needages to surface. So a compromise should be achieved by experienceand looking at how fast the air bubbles are surfacing.153Appendix B. Tissue Mimicking Materials• The temperature at which the molten plastic is kept to allow for theair bubbles to surface depends on the ratio of hardener or softener toplastic. Generally a harder plastic requires a higher temperature tobe 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 avery soft plastic and 180C for a very hard one.• When the liquid touches the cold metal surfaces of a mold, it willimmediately cool down and coagulate. The best practice would be toheat up the mold to the same temperature as liquid before pouringin the molten plastic. Otherwise, care should be taken in pouring themolten plastic. Too slow a rate would create a laminated phantombecause of the cooling down of the previously poured material. On theother had too fast a rate would trap the air bubbles in the phantom.• The ratio of softener or hardener to plastic volume should not exceedone because the created phantom would be too soft or too hard andvery difficult to handle. 2/3 limit would be a more conservative optioncreating a large elasticity interval (about one decade).• During the heating process of the liquid plastic, at the gel point solidmasses start to form inside the solution which would impede the stir-ring process. If the magnetic stirrer is used a lower stirring rate couldbe used until the solid masses are melted and a uniform liquid isachieved. Generally speaking the stirring should be kept to a mini-mum as it sucks in air bubbles. However without stirring it would bedifficult to maintain an even temperature distribution in the liquid.154References[1] T. J. Hall, M. Bilgen, M. F. Insana, and T. A. Krouskop, “Phantommaterials for elastography,” IEEE Transactions on Ultrasonics, Ferro-electrics, and Frequency Control, vol. 44, no. 6, pp. 1355–1365, 1997.[2] E. Madsen, G. Frank, T. Krouskop, T. Varghese, F. Kallel, andJ. Ophir, “Tissue-mimicking oil-in-gelatin dispersions for use in het-erogeneous 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 dis-persions 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 hetero-geneous elastography phantoms,” Physics in Medicine and Biology,vol. 50, p. 55975618, 2005.[5] C. M. Hassan and N. A. Peppas, “Structure and applications ofpoly(vinyl alcohol) hydrogels produced by conventional crosslinking orby 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, “Anthropomor-phic breast phantoms for testing elastography systems,” Ultrasound inMedicine and Biology, vol. 32, no. 6, pp. 857–874, 2006.[7] H. Asaba, E. Ohdaira, N. Masuzawa, and M. Ide, “Fundumental studyto develop bone-mimicking phantom,” Japanese Journal of AppliedPhysics, vol. 38, no. 5B, pp. 3412–3413, May 1999.[8] T. Hirai, E. Ohdira, N. Masuzawa, K. Itoh, and T. Matozaki, “Fun-damental study of bone-mimicking phantom using apetite,” JapaneseJournal of Applied Physics, vol. 40, no. 5B, pp. 3905–3906, May 2001.155Appendix B. Tissue Mimicking Materials[9] M. Belohlavek, V. B. Bartleson, and M. E. Zobitz, “Real-time strainrate imaging:validation of peak compression and expansion rates by atissue-mimicking phantom,” Echocardiography: A Jrnl. of CV Ultra-sound 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 testphantoms,” Ultrasound in Med. & Biol., vol. 29, no. 7, pp. 1053–1060,2003.[11] E. Brusseau, P. Delachartre, and D. Vray, “Axial strain imaging of ves-sel mimicking cryogel phantoms,” in Proceedings of Ultrasonics Sympo-sium 2000 IEEE, vol. 2. IEEE, October 2000, pp. 1817–1820.[12] K. C. Chu and B. K. Rutt, “Polyvinyl alcohol cryogel: An ideal phan-tom material for mr studies of arterial flow and elasticity,” MagneticResonance 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 char-acterization of a tissue-equivalent phantom for optical elastography,”Journal of Biomedical Optics, vol. 10, no. 4, p. 044020, July-August2005.[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 astissue phantoms in photoacoustic mammography,” Physics in Medicineand Biology, vol. 48, pp. 357–370, 2003.[15] K.Surry, H.Austin, A.Fenster, andPeters, “Poly(vinylalcohol)cryogelphantoms for use in ultrasound and mr imaging,” Physics in Medicineand Biology, vol. 49, p. 55295546, 2004.[16] W. D. D’Souza, E. L. Madsen, O. Unal, K. K. Vigen, G. R. Frank, andB. R. Thomadsen, “Tissue mimicking materials for a multi-imagingmodality 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 Con-trol, vol. 51, no. 24, pp. 410–419, April 2004.156Appendix B. Tissue Mimicking Materials[18] S. Girnyk, A. Barannik, E. Barannik, V. Tovostiak, A. Marusenko, andV. Volokhov, “The estimation of elasticity and viscosity of soft tissuesin vitro using the data of remote acoustic palpation,” Ultrasound inMed. and Biol., vol. 32, no. 2, p. 211219, 2006.[19] S. Inglis, K. V. Ramnarine, J. N. Plevris, and W. N. McDicken, “Ananthropomorphic tissue-mimicking phantom of the oesophagus for en-doscopic 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 ad-justable elastic and echographic properties,” in Proceedings of 2004IEEE International Ultrasonics, Ferroelectrics, and Frequency ControlJoint 50th Anniversary Conference. IEEE, 2004, pp. 1502–1505.[21] T. Kondo and H. Fujimoto, “Ultrasound tissue-mimicking materials us-ing oil gel and measurement of their characteristics,” Japanese Jounralof 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 performanceof 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, andD. Bryk, “Acoustic properties of polyvinyl chloride gelatin for use inultrasonography,” Radiology, pp. 215–216, 1984.[24] E. J. Chen, J. Novakofski, K. Jenkins, and W. D. O. Jr., “Ultrasoundelasticity measurements of beef muscle,” in Proc. of IEEE 1994 Ultra-sonics Symposium, 1994, pp. 1459–1462.[25] E. J. Chen, J. Novakofski, W. K. Jenkins, and W. D. O. Jr., “Youngsmodulus measurements of soft tissues with application to elasticityimaging,” IEEE Transcations on Ultrasonics, Ferroelectrics, and Fre-quency 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 tipgeometry of surgical needles: An assessment of force and deflection,”in Proc. of the 3rd European Medical and Biological Engineering Con-ference, Prague, Czech Republic, Nov. 2005.157Appendix CPhasor CalculationsConsider a scenario in which the tissue is continuously excited by a vibratingmechanical exciter. Practical limitations such as the sampling frequency,the finite time of observation, and the synchronization of the imaging andexcitation systems impose constraints on the (time) shape of the excitation.The excitation is usually chosen to be a periodic signal. Without a lossof generality, assume that the excitation is sinusoidal at a frequency of fe.When multiple frequencies exist in the excitation, the exact same analysisapplies for each of the components if linearity is assumed.In the steady state, the displacement of the point ¯x = (x,y,z) in thetissue at time t can be written as,ui(¯x,t) = Ai(¯x)exp(j2pifet)exp(jφi(¯x))+Ai(¯x)exp(−j2pifet)exp(−jφi(¯x)) (C.1)i = x,y,zwhere ¯u = (ux,uy,uz) is the displacement vector. Since the two terms onthe right hand side are complex conjugates, it is sufficient to use one for theanalysis. We call ¯U = (Ux,Uy,Uz) defined byUi(¯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 computeit, however, the displacements must be observed over an interval of time.Denoting the Fourier transform in the time variable by Ft{·}, ¯U can becomputed from ¯u,¯U(¯x) = Ft{¯u(¯x,t)}|ω=2pife (C.3)The phasor profile is used in the inversion techniques to estimate thetissue elasticity. A group of the techniques are based on the direct inversionof the wave equation. In some of these techniques, it is assumed that a com-ponent of the displacements satisfies the shear wave equation. For instance158Appendix C. Phasor Calculationsfor the y-component, one has∂∂t2uy(¯x,t) = c2s(¯x)∇2uy(¯x,t) (C.4)where cs(¯x) = radicalbigµ(¯x)/ρ is the classically known shear wave speed and ∇2is the Laplacian,∇2 = ∂2∂x2 +∂2∂y2 +∂2∂z2 (C.5)The equation for the phasor profile can be found by taking the Fouriertransform of (C.4) in the time variable, and evaluating it at the frequencyof excitation as in (C.3),−(2pife)2Uy(¯x) = cs(¯x)2∇2Uy(¯x) (C.6)The shear modulus can be computed fromµ(¯x) = ρc2s(¯x)= −ρ(2pife)2 Uy(¯x)∇2Uy(¯x)(C.7)Another technique which is used for the inversion is the phase gradienttechnique. Assume that a plane wave is propagating in the ˆn direction inthe medium. Assume that the coordinate axes have been chosen, so thatthe direction of particle motion for the wave is in the y-direction. They-component of the displacement can be written asuy(¯x,t) = Aexp(j2pi(fet− 1λˆn· ¯x))= Aexp(j2pi(fet))exp(−jkˆn· ¯x) (C.8)where λ is the wave length and k is the wave number,k = 2piλ = 2pifecs(C.9)Note that for brevity the conjugate term present in (C.2) has been droppedin (C.8). The phasor profile for this traveling wave becomes,Uy(¯x) = Aexp(−jkˆn· ¯x) (C.10)159Appendix C. Phasor CalculationsThe 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 numberk is first found. This is done by taking the gradient of the phase in the ˆndirection,k = −ˆn·∇∠Uy(¯x) (C.12)The shear wave speed is then given by,cs = − 2pifeˆn·∇∠Uy(¯x)(C.13)and the shear modulus by,µ = ρparenleftbigg 2pifeˆn·∇∠Uy(¯x)parenrightbigg2(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-12-31
pdf
Page Metadata
Item Metadata
Title | A wave equation approach to ultrasound elastography |
Creator |
Baghani, Ali |
Publisher | University of British Columbia |
Date | 2010 |
Date Issued | 2010-01-12T15:20:12Z |
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 |
Collection |
Electronic Theses and Dissertations (ETDs) 2008+ |
Date Available | 2010-01-12 |
Provider | Vancouver : University of British Columbia Library |
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 |
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