On the Identification of Mechanical Properties of Viscoelastic Materials by Hani Eskandari B.A.Sc., Sharif University of Technology, 2001 M.A.Sc., Sharif University of Technology, 2003 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University of British Columbia (Vancouver) March, 2009 c© Hani Eskandari 2009 Abstract Commonly used medical imaging techniques can render many properties of the anatomy or function, but are still limited in their ability to remotely measure tissue mechanical properties such as elasticity and viscosity. A remote and objective palpation function would help physicians in locating possible tumors or malignancies. The branch of medical imaging that char- acterizes tissues mechanical properties in a non-invasive manner has enjoyed increasing interest in the past two decades. The basic principle is to apply an excitation, such as tissue compression, to a region of interest and mea- sure the resulting tissue deformation. Tissue mechanical properties can then be inferred from the observed deformation at multiple locations in the re- gion, and the properties can be displayed as an image. If the excitation is dynamic, the deformation is considered as a motion field that varies in time and location over the region of interest. Ultrasound is particularly well suited for measuring motion fields due to its ability to image in real-time, low cost, low risk and ease-of-accessibility. The focus of this thesis is the estima- tion of the viscoelastic parameters such as Young’s modulus, viscosity and relaxation-time. For this purpose, a motion estimation method is proposed to measure axial tissue displacements from the peak of the ultrasound radio frequency signals. The displacements can be further processed to identify the mechanical properties. Two methods were developed: the first one is based on a one dimensional Voigt’s model of soft tissue and the second one is based on a finite element model. In the first method, a single frequency or wide-band excitation is applied to the tissue and the local relaxation-time is recovered from the phase difference between the strains or displacements. In this method, the elasticity can also be reconstructed from the magnitudes of the spectra. In the second approach, a novel dynamic finite element model is proposed for the incompressible soft materials. An inverse problem of viscoelasticity is solved iteratively to reconstruct the viscosity and elasticity based on a two or three dimensional model. The theoretical aspect of com- pressional elastography and longitudinal wave propagation is investigated. It is shown to be feasible to apply dynamic or transient compressional exci- tation to recover the dynamic properties of soft tissue. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . x Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Strain Imaging . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Static Elastography Techniques . . . . . . . . . . . . 4 1.2.3 Dynamic Elastography Techniques . . . . . . . . . . . 5 1.3 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . 8 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2 Tissue Strain Imaging Using a Wavelet Transform Based Peak Search Algorithm . . . . . . . . . . . . . . . . . . . . . . 17 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.1 Peak Detection Procedure . . . . . . . . . . . . . . . 21 2.2.2 Peak Matching . . . . . . . . . . . . . . . . . . . . . . 24 2.2.3 Strain Estimation . . . . . . . . . . . . . . . . . . . . 26 2.3 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 iii Table of Contents 2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.4.1 Numerical Results . . . . . . . . . . . . . . . . . . . . 27 2.4.2 Experimental Results . . . . . . . . . . . . . . . . . . 34 2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3 Viscoelastic Parameter Estimation Based on Spectral Analy- sis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2 Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.2.1 Linear Viscoelastic Model . . . . . . . . . . . . . . . . 51 3.2.2 Spectral Analysis and the Transfer Functions . . . . . 54 3.2.3 Power-Law Effect . . . . . . . . . . . . . . . . . . . . 54 3.3 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.3.1 Parameter Estimation . . . . . . . . . . . . . . . . . . 57 3.4 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.5.1 Actuator and Rheometer Design . . . . . . . . . . . . 61 3.5.2 Phantom Construction and Rheometry . . . . . . . . 64 3.5.3 Data Analysis and Results . . . . . . . . . . . . . . . 69 3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4 Viscoelastic Characterization of Soft Tissue from Dynamic Finite Element Models . . . . . . . . . . . . . . . . . . . . . . 82 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.2.1 Linear Viscoelastic Model . . . . . . . . . . . . . . . . 85 4.2.2 Finite Element Model of the Deformation . . . . . . . 86 4.3 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.3.1 Reconstruction Method . . . . . . . . . . . . . . . . . 87 4.3.2 Calculation of the Jacobian and the Hessian . . . . . 88 4.3.3 Practical Considerations . . . . . . . . . . . . . . . . 89 4.4 Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4.1 Two Dimensional Modeling . . . . . . . . . . . . . . . 90 4.4.2 Numerical Simulations . . . . . . . . . . . . . . . . . 94 iv Table of Contents 4.4.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . 102 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5 Model Validation for Longitudinal Wave Propagation in a Viscoelastic Environment . . . . . . . . . . . . . . . . . . . . . 117 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.2 Governing Equations in a Linear Viscoelastic Medium . . . . 118 5.3 Wave Equation and the Propagation Speed . . . . . . . . . . 121 5.3.1 Longitudinal Wave in a Finite Medium Approximat- ing a Thin Rod . . . . . . . . . . . . . . . . . . . . . 121 5.3.2 Longitudinal Wave Induced by a Point Source Exciter 123 5.4 Determining the Wave Speed from the Spectra of the Dis- placements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 5.4.1 Case (i) . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.4.2 Case (ii) . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.5.1 Dynamic Finite Element Analysis . . . . . . . . . . . 131 5.5.2 Simulation of the Wave Propagation . . . . . . . . . . 133 5.5.3 Model Comparison . . . . . . . . . . . . . . . . . . . 139 5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 6 Conclusions and Future Research . . . . . . . . . . . . . . . . 153 6.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.1.1 Strain Imaging . . . . . . . . . . . . . . . . . . . . . . 153 6.1.2 Spectral Viscoelastic Parameter Estimation . . . . . . 155 6.1.3 Finite Element Simulations . . . . . . . . . . . . . . . 156 6.1.4 Finite Element Inverse Problems . . . . . . . . . . . . 156 6.1.5 Justifying the Feasibility of Compressional Elastogra- phy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 6.1.6 Analysis of the Longitudinal Wave Propagation . . . 158 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 v Table of Contents Appendices A Wavelet Transforms and Peak Detection . . . . . . . . . . . 162 B Regularity of the Wavelet Transform . . . . . . . . . . . . . 165 C Transfer Function Parameterization . . . . . . . . . . . . . . 166 D Elasticity Estimation from the Transfer Functions . . . . . 168 E Simulating a Mass-Spring-Damper Network . . . . . . . . . 169 vi List of Tables 3.1 Rheometer Specifications . . . . . . . . . . . . . . . . . . . . 74 3.2 Power-law fit of relaxation-time . . . . . . . . . . . . . . . . . 74 5.1 Wave speeds for different boundary conditions . . . . . . . . . 137 vii List of Figures 1.1 A typical system for ultrasound elastography . . . . . . . . . 3 2.1 Typical RF A-line before and after compression . . . . . . . . 20 2.2 Wavelet transform modulus maxima and zero-crossings of the signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3 Block diagram of the peak matching procedure. . . . . . . . . 25 2.4 Estimated vs. real strain under different compressions for PSA and CC . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.5 Strain filter for PSA . . . . . . . . . . . . . . . . . . . . . . . 30 2.6 RMS errors for PSA . . . . . . . . . . . . . . . . . . . . . . . 32 2.7 Sensitivity for PSA . . . . . . . . . . . . . . . . . . . . . . . . 33 2.8 Effect of the sonographic noise on the PSA . . . . . . . . . . 35 2.9 PSA strain imaging of a simulated phantom with an inclusion 36 2.10 PSA strain imaging of a three-layered tissue mimicking phantom 38 2.11 PSA strain imaging of a tissue mimicking phantom with a hard inclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 A lumped mass-spring-damper model of soft tissue . . . . . . 50 3.2 The effect of the viscoelastic parameters on the transfer func- tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.3 The effect of the displacement SNR on the estimation of elas- ticity and relaxation-time . . . . . . . . . . . . . . . . . . . . 60 3.4 The effect of the LSQ kernel size on the estimation of elasticity and relaxation-time . . . . . . . . . . . . . . . . . . . . . . . . 62 3.5 The schematic of the experimental setup . . . . . . . . . . . . 65 3.6 Rheometry of the gelatin samples and the PVA sponge . . . . 66 3.7 Transfer function-based viscoelastic parameter estimation in a single line for a gelatin phantom with a viscous inclusion . . 68 3.8 Viscoelastic parameter estimation in a phantom with an elas- ticity and a viscosity inclusion using the transfer functions . . 71 4.1 A 2D finite elements mesh with rectangular elements . . . . . 91 viii List of Figures 4.2 Block diagram of the optimization procedure. . . . . . . . . . 92 4.3 The axial strain distribution around an elasticity and a vis- cosity inclusion due to different loading profiles. . . . . . . . . 93 4.4 The effect of the displacement SNR on the estimated para- meters using FEA . . . . . . . . . . . . . . . . . . . . . . . . 96 4.5 Diagram of a simulated phantom and a region-of-interest to test the effect of unknown boundary conditions . . . . . . . . 98 4.6 The effect of unknown boundary conditions in the reconstruc- tion of elasticity and viscosity using FEA. . . . . . . . . . . . 99 4.7 The effect of 3D deformations on finite element-based para- meter estimation when a 2D model is used for the inverse problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.8 The effect of uncertainty in the viscosity finite element model on the reconstruction of the viscosity and elasticity. . . . . . . 109 4.9 Finite element based viscosity and elasticity reconstruction of a gelatin phantom with an inclusion . . . . . . . . . . . . . . 110 5.1 Cylindrical coordinates . . . . . . . . . . . . . . . . . . . . . . 125 5.2 Displacement magnitude of the wave equation solution versus depth and frequency . . . . . . . . . . . . . . . . . . . . . . . 129 5.3 3D region geometry . . . . . . . . . . . . . . . . . . . . . . . . 132 5.4 Transfer functions for different boundary conditions . . . . . 135 5.5 Effect of the viscosity on the transfer function . . . . . . . . . 137 5.6 Transient step response for different boundary conditions . . 138 5.7 Illustration of longitudinal wave propagation for different bound- ary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 5.8 Comparison of the transfer functions of 1D, 2D and 3D models141 5.9 Comparison of the displacements and strains of 1D, 2D and 3D models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.10 Fitting different models to an experimental result . . . . . . . 145 A.1 Spectral characteristics of typical RF signal and the wavelet filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 ix List of Abbreviations 1D One Dimentional 2D Two Dimensional 3D Three Dimensional ARFI Acoustic Radiation Force Impulse CNR Contrast to Noise Ratio CC Cross Correlation CR Compression Ratio CT Computed Tomography ECG Electrocardiogram FEM Finite Element Method IR Infra Red KAVE Kinetic Acoustic Vitreo-retinal Examination LED Light Emitting Diode LSQ Least Squares MRE Magnetic Resonance Elastography MRI Magnetic Resonance Imaging MSD Mass-Spring-Damper PSA Peak Search Algorithm PSD Position Sensing Device PVA Polyvinyl Alcohol PVC Polyvinyl Chloride RF Radio Frequency RMS root mean square ROI Region of Interest RT Relaxation Time SNR Signal to Noise Ratio SSI Supersonic Shear Imaging STD Standard Deviation SWEI Shear Wave Elasticity Imaging TF Transfer Function WT Wavelet Transform x Acknowledgments I would like to thank Professor Salcudean and Professor Rohling for their valuable supervision and insight into my research during the course of this thesis. My special thanks to my parents whose unconditional support and love helped me accomplish this and many other milestones throughout my life. Many thanks and love to my wife, Zahra, the most beautiful result of my Ph.D., who has been with me in every step ever since we met half way through my research. xi Statement of Co-Authorship This thesis is based on several manuscripts, resulting from collaboration between multiple researchers. A version of Chapter 2 appeared in the IEEE Transactions on Ultra- sonics, Ferroelectrics and Frequency Control in June, 2007. This paper was co-authored with Dr. Tim Salcudean and Dr. Robert Rohling. The author’s contribution in that paper was developing the idea, numerical simulation, phantom fabrication, experimental evaluation of the algorithms and writ- ing the manuscript. Dr. Salcudean and Dr. Rohling assisted with their suggestions and editing the manuscript. A version of Chapter 3 appeared in the IEEE Transactions on Ultra- sonics, Ferroelectrics and Frequency Control in July, 2008. This paper was co-authored with Dr. Tim Salcudean and Dr. Robert Rohling. The au- thor’s contribution in that paper was developing the equations culminating in the relaxation-time estimation, numerical simulation, phantom fabrica- tion, experimental evaluation of the algorithms, designing the position sen- sor, calibration of the experimental set-up, software development to control the rheometer and writing the manuscript. Dr. Salcudean designed the rheometer platform and suggested the author to investigate the phase and other aspects of the transfer functions that can potentially relate to vis- coelastic parameters. Dr. Rohling and Dr. Salcudean helped with their numerous suggestions in the course of enhancing the data acquisition and parameter identification procedures as well as editing the manuscript. A version of Chapter 4 appeared in the Physics in Medicine and Biology in November, 2008. This paper was co-authored with Dr. Tim Salcudean, Dr. Robert Rohling and Dr. Jacques Ohayon. The author’s contribution in that paper was formulating the inverse problem and development of the solution strategy, developing the damping model for the finite element analy- sis, numerical simulation, phantom fabrication, experimental evaluation of the algorithms, designing the position sensor, calibration of the experimen- tal set-up, software development to control the rheometer and writing the manuscript. Dr. Salcudean suggested the finite element method as a tool to investigate the inverse problem of viscoelasticity. Dr. Salcudean and xii Statement of Co-Authorship Dr. Rohling assisted through the course of concept development and helped in editing the manuscript. Dr. Ohayon suggested and showed the way to developed the new definition of the damping matrix for incompressible vis- coelastic materials. A version of Chapter 5 has been submitted to Physics in Medicine and Biology in March, 2009. This paper is co-authored with Ali Baghani, Dr. Tim Salcudean and Dr. Robert Rohling. The author’s contribution in that paper was expanding the wave equation and analyzing the effect of differ- ent boundary conditions as an original contribution, numerical simulations, experimental validation of the theories and writing the manuscript. Ali Baghani helped with developing the concept and suggested exploiting the wave equation for understanding the propagation of the longitudinal waves as well as editing the manuscript. Dr. Salcudean and Dr. Rohling helped with editing the manuscript and general research directions. xiii Chapter 1 Introduction 1.1 Motivation Seeing through the human body in order to delineate changes and patholo- gies has been a challenge for the modern science. Several techniques have been developed to depict different tissue parameters [1]. The foremost is the X-ray imaging which projects the volumetric X-ray absorption coefficients on a two dimensional image. Three dimensional mapping of these absorption factors has been developed on the first hand in the late 70s which was the first step toward development of the computer tomography (CT) technology. Later on and especially in the last decade, with the introduction of spiral and multi-slice scanners, promising improvements regarding the speed and the diagnostic aspects have been introduced in the field of the computed tomog- raphy [2]. Magnetic Resonance Imaging (MRI) is also another powerful tool that reveals the proton density together with the damping time-constant of the dipoles’ spin inside tissue. The least expensive and one of the handiest and least harmful techniques of imaging the human body is ultrasound. Any acoustic or ultrasonic beam passing through a medium will be partially re- flected on the boundaries between tissues, depending on the difference of the acoustic impedances between the two regions. Using the ultrasound tech- nique, one can find the boundaries and scatterers that exist inside a region of interest [3]. The motivation for establishing a new modality of medical imaging, be- sides the cost and slowness of some of the above imaging techniques is that in the presence of many cancerous tumors and malignant carcinoma, con- ventional methods are sometimes unable to discriminate between the patho- logical and healthy tissues while an expert physician can still localize them by simple fingertip palpation. In order to illustrate the stiffness of tissues, methods have been devised to compute the Young’s Modulus or the strain distribution inside an imaging window. These methods are based on taking images of the medium under normal and compressed states and finding the displacements of different points or patches and finally analyzing them to extract the desired parameters. Generally, this technique is being referred 1 Chapter 1. Introduction to as elastography in the literature [4]. 1.2 Background Different imaging modalities can be used in order to estimate the tissue deformation, among which ultrasound and MRI are the most popular ones. On the one hand, ultrasound, due to its low cost, portability and relatively high frame-rate, is the most widely used modality, which has made possible the development of several elastography techniques [4–16]. On the other hand, in magnetic resonance elastography (MRE) it is possible to obtain a 3D map of the tissue deformation with a higher signal-to-noise ratio and overall resolution than ultrasound [17–21]. A typical system to estimate the mechanical properties of soft tissue with ultrasound is depicted in Fig. 1.1. 1.2.1 Strain Imaging When a medium undergoes compression or shearing, the interior region will experience a nonuniform deformation profile. These displacements can be detected by several algorithms of motion detection and the strain is ex- tracted by taking the gradient or local slope of the displacements. Because of the high axial resolution of ultrasound radio frequency (RF) signals, more accurate estimates can be obtained from processing RF data compared to the image processing techniques based on ultrasound B-mode images. Generally, in the static case, a pre-compression and a post-compression RF signal are compared to extract the 1D motion in that line. In the dynamic case where excitation is continuously applied to the tissue, each RF frame can be compared to the previous frame to estimate relative motions and then integrate. This way, due to non-linearities in the 1D projected model of tissue deformation, the overall displacements may be subject to a high amount of drift. In another approach, absolute estimation is possible through using one frame as the reference frame and comparing all other frames to that one. While the drift will be reduced this way, the accuracy of the estimations will be compromised as the similarity between pre- and post-compression signals decreases as imaging time increases. Several methods are proposed in the literature, base on mutual features or correlation between the two signals. Time domain cross correlation is a conventional method that divides the two signals into equally spaced, over- lapping blocks and correlates each block in one signal with the block in the other signal to find the motion associated with that block [22–24]. Consider- ing that the transformation between the pre- and post-compression signals is 2 Chapter 1. Introduction t= t1 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t = t0 Parameter Estimation Motion Estimation Ultrasound Images Tissue vibration or Compression Displacements (u) Gradient Mechanical Parameters Strains Boundary Conditions Figure 1.1: A typical system for ultrasound elastography of a tissue- mimicking phantom. The excitation can be static or dynamic. The B-mode images or RF data can be used to measure the deformations. Strains can be calculated from the displacements using gradient or least-square approaches. Depending on the model, some information of the boundary conditions may be required to extract the mechanical parameters. 3 Chapter 1. Introduction not limited to block motions, but also overall stretching is needed to recover the original signals, accuracy can be improved by global stretching of the compressed signal prior to motion estimation [25]. Adaptive stretching can also be incorporated with a cross-correlation algorithm to add another de- gree of freedom for matching corresponding blocks [26]. The higher accuracy of adaptive stretching, however, comes at the expense of its low computa- tional efficiency. Since the time delay between blocks shows itself as a phase lag in the frequency domain, finding this lag has been another way to esti- mate the motion. This method, known as zero-phase seeking algorithm, is fast but less accurate in small compressions [27, 28]. Viola et. al presented a spline-based motion estimation technique that directly determines sub- sample time delay estimates from the RF data [29,30]. The algorithm uses cubic splines to obtain a continuous representation of a reference signal. An analytical representation of a matching function is then computed between the reference and a compressed signal. The matching function is optimized to obtain the time delays. Low lateral resolution of ultrasound limits the achievable accuracy for es- timating the lateral displacements [31,32]. The available techniques mostly rely on lateral or two-dimensional (2D) feature matching or correlation be- tween sequential RF frames. Three dimensional displacement estimation by means of ultrasound is possible through translation, rotation or fanning of the 2D transducer array. In [33] and [34] a feature tracking method was investigated which can be suitable for 3D blood flow velocity measurement. Given the low frame-rate and very limited elevational resolution of the 3D ultrasound systems, MRI is the preferred technology to estimate the 3D displacement profile in soft tissues [35]. 1.2.2 Static Elastography Techniques In the 3D static case, the partial differential equations of the continuum that describe the relationship between the stress and strain can be analyzed in 3D to solve for the anisotropic parameters of the medium. Assuming a nearly incompressible and isotropic medium significantly reduces the number of elasticity parameters and leaves the Young’s modulus as the only unknown. If the problem is discretized and formulated such that the Young’s moduli are the unknowns, with known displacements and strains, a forward solution technique can be applied to estimate the parameters [36–41]. Most of the methods in this category require the 2D displacement and strain fields as well as their first and second derivatives in the region of interest. Solving the elasticity distribution from the 2D or 3D strain distribution is a challeng- 4 Chapter 1. Introduction ing inverse problem in the field of elastography. A given strain distribution in a 3D medium can be a result of infinitely many elasticity distribution cases [42]. Therefore, additional information such as boundary conditions and assumptions on tissue mechanical properties is necessary to obtain a unique solution. For example, most of the direct elasticity reconstruction methods are limited to 2D, isotropic, elastic, incompressible media. Sumi et. al proposed a method to estimate the shear modulus distribution for such a medium [37]. In that method, a linear set of equations were derived with spatial derivatives of the shear modulus as the unknowns and the 2D strain data as the coefficients of the equations. These equations were solved ana- lytically to obtain the spatial derivatives of the shear modulus distribution in the region of interest (ROI). The elasticity distribution was then found by integration. Skovoroda et. al used a two-step procedure to solve for the elas- ticity distribution in the presence of an inclusion with sharp boundaries [38]. First, the stress continuity condition was used to find the discontinuities in the shear modulus distribution. The discontinuities were then used as the regional boundaries for solving the partial differential equations of the con- tinuum. Due to the low lateral resolution of ultrasound, the estimated lateral displacements have low signal-to-noise ratios (SNR). Once the low SNR of the lateral displacements has been overcome and the boundary forces have been measured, a reliable elasticity distribution may be obtained [41]. In another approach, the inverse problem of elasticity is solved in the sense that a specific functional is minimized [43–49]. Usually, this func- tional is considered to be the quadratic norm of the difference between the measured displacements and the displacements resulting from an assumed distribution of elasticity. The medium can be discretized using finite differ- ences or finite element method (FEM), and iterative strategies can be used to solve for the unknown elastic parameters. The exact analytical solu- tions for the inverse problem of elasticity and other quasi-static parameters have been investigated in [50], where it was shown that the problem of de- termining the shear modulus becomes unstable for nearly incompressible materials. It was shown that the uniqueness of the solution is highly de- pendent on the regularization, information about the boundary conditions, prior knowledge of the elasticity on the boundaries and the incompressibility of the medium [42,51,52]. 1.2.3 Dynamic Elastography Techniques Applying dynamic excitation to the tissue can reveal additional properties of the medium such as wave speed, viscosity and steady-state harmonic be- 5 Chapter 1. Introduction havior. Conventional compressional elastography involves applying a static compression to the tissue externally and solving for the Young’s modulus or interpreting the strain images. However, applying a transient excita- tion in the same manner, known as vibro-elastography, can potentially yield viscoelastic properties of the medium [15, 53]. The excitation in this case can be a single-frequency or a wide-band vibration to the surface of the tis- sue. Different methods of identification can be utilized, including parametric modeling, non-parametric spectral analysis, FEM-based inverse problems and analysis of the wave front. The damping, which if often due to the viscous effects in soft tissue, results in a time-constant decay in the stress relaxation and creep tests on tissue-mimicking phantoms. Insana et. al pro- posed a method to recover the time-constants from the stress relaxation and creep response of soft tissues imaged by ultrasound [11,12]. In sonoelasticity imaging a low-frequency vibration (around 100Hz) is applied to tissue, while being imaged with Doppler ultrasound [13]. A small stiff inhomogeneity in the surrounding tissue appears as a disturbance in the normal vibration eigen-mode pattern and can be detected thereby. Alternatively, a shear excitation can be applied to the tissue and the re- sulting deformation can be analyzed. The shear wave can be induced either by an external surface exciter or through acoustic radiation force impulse (ARFI) which uses high frequency ultrasound waves to create a force field inside the medium [54]. As a result, shear wave starts to propagate perpen- dicular to the direction of excitation. Compared to the speed of purely lon- gitudinal waves, such as ultrasound that travels at approximately 1540m/s in soft tissue, the shear wave speed is generally less than 10m/s. Therefore, it is feasible to follow the shear wave front with the help of ultrasound or MRI. Shear wave elasticity imaging (SWEI) [5] and ARFI imaging [6] are two techniques that are based on focusing an ultrasound beam for a short period of time on a location inside tissue and monitoring the resulting dis- placements at the focus point. In vibro-acoustography two ultrasound beams of slightly different frequencies are emitted to the tissue which generate si- nusoidal vibrations in a portion of the object [7, 8]. The displacements, being a function of the viscoelastic properties of the tissue, are measured by a hydrophone or another ultrasound transducer. The dispersion produces low-frequency shear wave for which the speed and attenuation can be asso- ciated with the shear elasticity and viscosity of the medium. In supersonic shear imaging (SSI) the point of focus of the ARFI is moved at supersonic speed and the resulting shear wave is ultrasonically tracked [9]. In kinetic acoustic vitreo-retinal examination (KAVE) the radiation force is used to generate a steady state stress field within the body of the vitreous of the eye 6 Chapter 1. Introduction and the ultrasonically imaged displacements are analyzed to characterize the viscoelastic properties [10]. 1.3 Thesis Objectives In this thesis new algorithms for imaging the viscoelastic properties of soft tissue are developed. The hypothesis that viscosity and elasticity can be estimated from the dynamic tissue response to a transient longitudinal or shear excitation will be investigated. To explore this hypothesis, the follow- ing objectives are defined in this thesis: 1. Accurate and efficient estimation of the axial motion as a result of tissue compression, using ultrasound radio-frequency data. 2. Robust estimation of the viscosity and elasticity from the time-dependent displacements, based on a 1D model of mass-spring-damper elements, by analyzing the transfer functions between local tissue displacements or strains. 3. Finite element modeling of soft tissue and viscoelastic parameter iden- tification based on a dynamic finite element model of the tissue. 4. Justifying the use of longitudinal excitation to image soft tissue defor- mation with ultrasound when proper boundary conditions are applied. 5. Constructing tissue-mimicking phantoms with controllable viscoelastic properties for validation purposes. In the course of achieving the primary objectives of this thesis, the fol- lowing contributions were made: • A strain imaging method was developed that uses the wavelet trans- form to estimate the axial motion of the medium from the location of peaks of the RF signals. • A method was developed to recover the relaxation-time of soft tissue in a vibro-elastography exam from the phase of the displacement or strain spectra based on a 1D model. • Nonlinear viscous effects in the form of a frequency-dependent power- law were observed at low-frequency in tissue-mimicking phantoms and a method was proposed to estimate the corresponding parameters from the spectral data. 7 Chapter 1. Introduction • Tissue simulation was performed with dynamic 2D and 3D finite ele- ment models. A new formulation of the damping matrix was proposed that is coherent with the Voigt’s viscoelastic model of soft tissue. • An inverse problem was formulated in a finite element framework to reconstruct the elasticity and viscosity of the medium. A method was proposed to estimate the parameters from harmonic displacements by minimizing a specific functional. • Gelatin-based tissue mimicking phantoms were constructed with dif- ferent embedded elasticity and viscosity inclusions in order to evaluate the performance of the estimation algorithms. The viscosity and elas- ticity of the inclusions could be independently controlled in such a way that the inclusion could have the same elasticity as the background but a different viscosity or vice versa. • The effect of the boundary conditions and excitation on the propaga- tion of longitudinal waves were analyzed. Through expansion of the wave equation for an incompressible viscoelastic medium, it was proved that the longitudinal wave can propagate at a speed of a few meters per second which makes them traceable by means of ultrasound. 1.4 Chapter Summary The overall goal of the research undertaken in this thesis was to develop al- gorithms for viscoelastic characterization of soft tissue and tissue-mimicking materials. The thesis presented here is written in the manuscript-based style, as permitted by the Faculty of Graduate Studies at the University of British Columbia. In the manuscript-based thesis, each chapter represents an indi- vidual work that has been published, submitted or prepared for submission to a peer reviewed publication. Each chapter is self-sustained in the sense that it includes an introduction to the work presented in that chapter, the methodology, results and discussion. The references are summarized in the bibliography found at the end of each chapter. The appendices pertaining to each chapter are presented at the end of the thesis. In Chapter 2 a new method is proposed to estimate the motion and the relative local compression between two successive ultrasound RF signals un- der different compression states. The algorithm uses the continuous wavelet transform to locate the peaks in the RF signals. The estimated peaks in the pre- and post-compression signals are assigned to each other by a peak 8 Chapter 1. Introduction matching technique with the goal of minimizing the number of false matches. The method allows local shifts of the tissue to be estimated. The method has been tested in one-dimensional simulations and phantom experiments. The signal-to-noise ratio and the rms error are shown to be better than for the standard cross-correlation method (CC). The new estimator remains unbiased for up to 10% strain which is a larger range than that of CC. The maximum signal-to-noise ratio is 3 times as high as that of the conventional CC method, showing higher sensitivity as well. The method is computation- ally efficient, achieving 0.7 msec/RF-line on a standard personal computer. In Chapter 3 a novel robust estimator for the relaxation-time of soft tissue is proposed. The main novelty is in the use of the phase of trans- fer functions calculated from a time series of strain measurements at mul- tiple locations. Computer simulations with simulated measurement noise demonstrate the feasibility of the approach. An experimental apparatus and software were developed to confirm the simulations. The setup can be used both as a rheometer to characterize the overall mechanical properties of a material or as a vibro-elastography imaging device using an ultrasound system. The algorithms were tested on tissue-mimicking phantoms specif- ically developed to exhibit contrast in elasticity and relaxation-time. The phantoms were constructed using a combination of gelatin and a polyvinyl alcohol sponge to produce the desired viscoelastic properties. The tissue parameters were estimated and the elasticity and relaxation-time of the ma- terials have been used as complementary features to distinguish different materials. The estimation results are consistent with the rheometry, veri- fying that the relaxation-time can be used as a complementary feature to elasticity to delineate the mechanical properties of the phantom. An iterative solution to the inverse problem of elasticity and viscosity is proposed in Chapter 4. A new dynamic finite element model that is consistent with known rheological models has been derived to account for the viscoelastic changes in soft tissue. The model assumes known lumped masses at the nodes, and comprises two vectors of elasticity and viscosity parameters that depend on the material elasticity and viscosity distribu- tions, respectively. Using this deformation model and the observed dynamic data for harmonic excitation, the inverse problem is solved to reconstruct the viscosity and elasticity in the medium by using a Gauss-Newton based approach. As in other inverse problems, previous knowledge of the parame- ters on the boundaries of the medium is necessary to assure uniqueness and convergence and to obtain an accurate map of the viscoelastic properties. The sensitivity of the solutions to noise, model and boundary conditions has been studied through numerical simulations. Experimental results are also 9 Chapter 1. Introduction presented. The viscosity and elasticity of a gelatin-based phantom with an inclusion of known properties have been reconstructed and have been shown to be close to the values obtained using standard rheometry. In Chapter 5 the effect of boundary conditions and excitation on the speed of the longitudinal wave in a finite medium are investigated. It is shown by analyzing the 3D wave equation and through finite element sim- ulations that a compressional excitation may travel at a very high speed or at speeds as low as that of the shear wave. As a result, besides the shear wave elastography, tracking the longitudinal waves in ultrasound elastogra- phy can be used as another method to investigate the viscoelastic behavior of the materials. A 3D FEM model is adopted as a reference for the analysis and validation of the theoretical results, while the error produced by sim- plified 2D FEM or even 1D models are explored. Transient and harmonic simulations are performed. A linear viscoelastic medium is assumed and the effect of viscosity on the displacement spectrum is discussed. In Chapter 6 the results of the collected works are related to one another and a unified goal of the thesis is discussed. The strengths and weaknesses of the research are then presented, along with future directions for research. 10 References [1] J. Beutel, H. Kundel and R.L. Van Metter, Handbook of Medical Imag- ing, Volume 1, SPIE Press, 2000. [2] M. Prokop and M. Galanski, Spiral and Multislice Computed Tomogra- phy of the Body, Georg Thieme Verlag, 2003. [3] D.A. Christensen, Ultrasonic Bioinstrumentation, Wiley, 1988. [4] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastog- raphy: a quantitative method for imaging the elasticity of biological tissues,” Ultrasonic Imaging, vol. 13, no. 2, pp. 111–34, April 1991. [5] A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov, “Shear wave elasticity imaging: a new ultrasonic tech- nology of medical diagnostics,” Ultrasound in Medicine and Biology, vol. 24, no. 9, pp. 1419–1435, 1998. [6] K. Nightingale, M. Palmeri, R. Nightingale, and G. Trahey, “On the feasibility of remote palpation using acoustic radiation force,” Journal of the Acoustical Society of America, vol. 110, no. 1, pp. 625–634, 2001. [7] M. Fatemi and J. F. Greenleaf, “Vibro-acoustography: An imaging modality based on ultrasound-stimulated acoustic emission,” Proceed- ings of the National Academy of Science, U.S.A., vol. 96, no. 12, pp. 6603–6608, June 1999. [8] M. Fatemi, A. Manduca, and J. F. Greenleaf, “Imaging elastic prop- erties of biological tissues by low-frequency harmonic vibration,” Pro- ceedings of the IEEE, vol. 91, no. 10, pp. 1503–19, Oct. 2003. [9] J. Bercoff, M. Tanter, and M. Fink, “Supersonic shear imaging: A new technique for soft tissue elasticity mapping,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 51, no. 4, pp. 396–409, 2004. 11 Chapter 1. Introduction [10] W. F. Walker, F. J. Fernandez, and L. A. Negron, “A method of imag- ing viscoelastic parameters with acoustic radiation force.” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1437–1447, Jun 2000. [11] M. F. Insana, C. Pellot-Barakat, M. Sridhar, and K. K. Lindfors, “Vis- coelastic imaging of breast tumor microenvironment with ultrasound,” Journal of Mammary Gland Biology and Neoplasia, vol. 9, no. 4, pp. 393–404, Oct. 2004. [12] M. Sridhar, J. Liu, and M. F. Insana, “Viscoelasticity imaging using ultrasound: parameters and error analysis.” Physics in Medicine and Biology, vol. 52, no. 9, pp. 2425–2443, May 2007. [13] L. Gao, K. J. Parker, S. K. Alam, D. Rubens and R. M. Lerner, “The- ory and application of sonoelasticity imaging” International Journal of Imaging Systems and Technology, vol. 8, no. 1, pp. 104–109 1997. [14] J. Ophir, I. Cespedes, B. Garra, H. Ponnekanti, Y. Huang, and N. Mak- lad, “Elastography: ultrasonic imaging of tissue strain and elastic mod- ulus in vivo,” European Journal of Ultrasound, vol. 3, no. 1, pp. 49–70, Jan. 1996. [15] H. Eskandari, S. Salcudean, and R. Rohling, “Viscoelastic parameter estimation based on spectral analysis,” IEEE Transactions on Ultrason- ics, Ferroelectrics and Frequency Control, vol. 55, no. 7, pp. 1611–1625, July 2008. [16] J. Ophir, S. K. Alam, B. Garra, F. Kallel, E. Konofagou, T. Korouskop, T. Varghese, “Elastography: ultrasonic estimation and imaging of the elastic properties of tissues”, Proc. Inst. Mech. Eng. [H], vol. 213, no. 3, pp. 203–233, May 2007. [17] R. Sinkus, M. Tanter, S. Catheline, J. Lorenzen, C. Kuhl, E. Son- dermann, and M. Fink, “Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography.” Magnetic Reso- nance in Medicine, vol. 53, no. 2, pp. 372–387, Feb 2005. [18] E. Ehman, P. Rossman, S. Kruse, A. Sahakian, and K. Glaser, “Vi- bration safety limits for magnetic resonance elastography,” Physics in Medicine and Biology, vol. 53, no. 4, pp. 925 – 935, 2008. [19] M. A. Green, L. E. Bilston, and R. Sinkus, “In vivo brain viscoelas- tic properties measured by magnetic resonance elastography,” NMR in Biomedicine, vol. 21, no. 7, pp. 755 – 764, 2008. 12 Chapter 1. Introduction [20] A. J. Romano, J. A. Bucaro, B. H. Houston, J. L. Kugel, P. J. Ross- man, R. C. Grimm, and R. L. Ehman, “On the feasibility of elastic wave visualization within polymeric solids using magnetic resonance elastography,” Journal of the Acoustical Society of America, vol. 116, no. 1, pp. 125 – 132, 2004. [21] I. Sack, B. Beierbach, U. Hamhaber, D. Klatt, and J. Braun, “Non- invasive measurement of brain viscoelasticity using magnetic resonance elastography,” NMR in Biomedicine, vol. 21, no. 3, pp. 265 – 271, 2008. [22] M. A. Lubinski, S. Y. Emelianov, and M. O’Donnell, “Speckle tracking methods for ultrasonic elasticity imaging using short-time correlation,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Con- trol, vol. 46, no. 1, pp. 82–96, Jan. 1999. [23] S. G. Foster, P. M. Embree, and W. D. J. O’Brien, “Flow velocity pro- file via time-domain correlation: Error analysis and computer simula- tion,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 37, no. 3, pp. 164–175, May 1990. [24] R. Zahiri-Azar and S. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, Oct. 2006. [25] T. Varghese and J. Ophir, “Enhancement of echo-signal correlation in elastography using temporal stretching.” IEEE Transactions on Ultra- sonics, Ferroelectrics and Frequency Control, vol. 44, no. 1, pp. 173–180, Jan. 1997. [26] S. K. Alam, J. Ophir, and E. E. Konofagou, “An adaptive strain estima- tor for elastography,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 45, no. 2, pp. 461–72, March 1998. [27] A. Pesavento, A. Lorenz, and H. Ermert, “Phase root seeking and the cramer-rao-lower bound for strain estimation,” in Proceedings of IEEE Ultrasonics Symposium, vol. 2, 1999, pp. 1669–72. [28] A. Pesavento, C. Perrey, M. Krueger, and H. Ermert, “A time-efficient and accurate strain estimation concept for ultrasonic elastography using iterative phase zero estimation,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 46, no. 5, pp. 1057–67, Sep. 1999. 13 Chapter 1. Introduction [29] F. Viola, M. A. Ellis, and W. F. Walker, “Time-domain optimized near-field estimator for ultrasound imaging: Initial development and results,” IEEE Transactions on Medical Imaging, vol. 27, no. 1, pp. 99–110, 2008. [30] F. Viola and W. F. Walker, “Computationally efficient spline-based time delay estimation,” IEEE transactions on ultrasonics, ferro- electrics, and frequency control, vol. 55, no. 9, pp. 2084–2091, 2008. [31] F. Kallel, T. Varghese, J. Ophir, and M. Bilgen, “Nonstationary strain filter in elastography: Part ii. lateral and elevational decorrelation,” Ultrasound in Medicine and Biology, vol. 23, no. 9, pp. 1357 – 1369, 1997. [32] R. Righetti, S. Srinivasan, and J. Ophir, “Lateral resolution in elastog- raphy,” Ultrasound in Medicine and Biology, vol. 29, no. 5, pp. 695 – 704, 2003. [33] G. R. Bashford and O. T. von Ramm, “Ultrasound three-dimensional velocity measurements by feature tracking,” IEEE Transactions on Ul- trasonics, Ferroelectrics, and Frequency Control, vol. 43, no. 3, pp. 376 – 384, 1996. [34] J. Kuo and O. T. Von Ramm, “Three-dimensional motion measure- ments using feature tracking,” IEEE Transactions on Ultrasonics, Fer- roelectrics, and Frequency Control, vol. 55, no. 4, pp. 800 – 810, 2008. [35] T. S. Denney and J. L. Prince, “Reconstruction of 3-D left ventricular motion from planar tagged cardiac MR images: an estimation theoretic approach,” IEEE Transactions on Medical Imaging, vol. 14, no. 4, pp. 625–635, Dec. 1995. [36] K. R. Raghavan and A. E. Yagle, “Forward and inverse problems in elasticity imaging of soft tissues,” IEEE Transactions on Nuclear Sci- ence, vol. 41, no. 4 pt 1, pp. 1639–1645, 1994. [37] C. Sumi, A. Suzuki, and K. Nakayama, “Estimation of shear modulus distribution in soft tissue from strain distribution,” IEEE Transactions on Biomedical Engineering, vol. 42, no. 2, pp. 193–202, Feb. 1995. [38] A. R. Skovoroda, S. Y. Emelianov, and M. O’Donnell, “Tissue elasticity reconstruction based on ultrasonic displacement and strain images,” IEEE transactions on ultrasonics, ferroelectrics, and frequency control, vol. 42, no. 4, pp. 747–765, 1995. 14 Chapter 1. Introduction [39] J. Bishop, A. Samani, J. Sciarretta, and D. B. Plewes, “Two- dimensional mr elastography with linear inversion reconstruction: methodology and noise analysis.” Phys Med Biol, vol. 45, no. 8, pp. 2081–2091, Aug 2000. [40] Y. Zhu, T. J. Hall, and J. Jiang, “A finite-element approach for young’s modulus reconstruction,” IEEE Transactions on Medical Imag- ing, vol. 22, no. 7, pp. 890–901, Sep. 2003. [41] E. Park and A. M. Maniatty, “Shear modulus reconstruction in dynamic elastography: time harmonic case.” Phys Med Biol, vol. 51, no. 15, pp. 3697–3721, Aug 2006. [42] P. E. Barbone and N. H. Gokhale, “Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions,” Inverse Problems, vol. 20, no. 1, pp. 283–96, Feb 2004. [43] F. Kallel and M. Bertrand, “Tissue elasticity reconstruction using linear perturbation method,” IEEE Transactions on Medical Imaging, vol. 15, no. 3, pp. 299–313, Jun 1996. [44] D. Fu, S. F. Levinson, S. M. Gracewski, and K. J. Parker, “Non-invasive quantitative reconstruction of tissue elasticity using an iterative forward approach,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1495– 509, Aug 2000. [45] M. I. Miga, “A new approach to elastography using mutual information and finite elements,” Physics in Medicine and Biology, vol. 48, no. 4, pp. 467–80, 02/21 2003. [46] 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, no. 13, pp. 2955–2974, Jul. 2004. [47] M. M. Doyley, E. E. V. Houten, J. B. Weaver, S. Poplack, L. Duncan, F. Kennedy, and K. D. Paulsen, “Shear modulus estimation using par- allelized partial volumetric reconstruction.” IEEE Trans Med Imaging, vol. 23, no. 11, pp. 1404–1416, Nov. 2004. [48] Y. Liu, L. Sun, and G. Wang, “Tomography-based 3-d anisotropic elas- tography using boundary measurements,” IEEE Transactions on Med- ical Imaging, vol. 24, no. 10, pp. 1323 – 1333, 2005. 15 Chapter 1. Introduction [49] A. S. Khalil, R. C. Chan, A. H. Chau, B. E. Bouma, and M. R. K. Mofrad, “Tissue elasticity estimation with optical coherence elastogra- phy: Toward mechanical characterization of in vivo soft tissue,” Annals of Biomedical Engineering, vol. 33, no. 11, pp. 1631–1639, 2005. [50] P. E. Barbone and A. A. Oberai, “Elastic modulus imaging: some exact solutions of the compressible elastography inverse problem.” Phys Med Biol, vol. 52, no. 6, pp. 1577–1593, Mar 2007. [51] P. E. Barbone and J. C. Bamber, “Quantitative elasticity imaging: what can and cannot be inferred from strain images.” Phys Med Biol, vol. 47, no. 12, pp. 2147–2164, Jun 2002. [52] J. R. McLaughlin and J.-R. Yoon, “Unique identifiability of elastic parameters from time-dependent interior displacement measurement,” Inverse Problems, vol. 20, no. 1, pp. 25–45, 2004. [53] E. Turgay, S. E. Salcudean, and R. Rohling, “Identifying the mechan- ical properties of tissue by ultrasound strain imaging,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 221–235, 2006. [54] O. V. Rudenko, A. P. Sarvazyan and S. Y. Emelianov, “Acoustic radi- ation force and streaming induced by focused nonlinear ultrasound in a dissipative medium”, J. Acoustical Society America, vol. 99, no. 5, pp. 2791-2798, May 1996. 16 Chapter 2 Tissue Strain Imaging Using a Wavelet Transform Based Peak Search Algorithm 1 2.1 Introduction Motion estimation in ultrasound images is one of the most important prob- lems in elastography. Estimating local displacements in single radio fre- quency (RF) lines is necessary to produce strain images by taking the gra- dient; it is also a necessary part of solving the inverse problem of elastic- ity [1–8]. Both strain estimation and inverse problems are very sensitive to input data errors. Tissue motion estimation from RF lines is also useful for reducing noise in ultrasound B-mode images [9]. Several methods of robust motion estimation have been introduced in the literature so far [10–18]. Maurice et al. utilized a Lagrangian speckle model to formulate the endovascular elastography problem in a polar coordinate system [10]. Sumi took advantage of the RF-echo phase information in 2-D to measure the displacement profile in a work-plane [11]. Pesavento et al. used a frequency-domain approach that measures the time delay between two signals as a phase difference. The phase difference is detected by taking the inverse Fourier Transform of the quotient of the two spectrums [12]. This method is very efficient computationally, but is of limited accuracy due to sampling. A more accurate version of this algorithm has also been proposed that iteratively considers the time shifts between the two signals to overcome problems such as aliasing [13]. The most widely used motion estimation method is the cross-correlation (CC) technique, in which the local relative compression between two RF lines 1A version of this chapter has been published. H. Eskandari, S. E. Salcudean, and R. Rohling, “Tissue Strain Imaging Using a Wavelet Transform Based Peak Search Algo- rithm”. IEEE Trans. on Ultras. Ferr. Freq. Cont., vol. 54, no. 6, pp. 1118-1130, Jun. 2007 17 Chapter 2. Peak Search Strain Imaging is determined by maximizing the cross-correlation between corresponding signal blocks [14, 15]. In this method, a compressed signal block is usu- ally selected by using a masking window. This block is swept over the uncompressed signal and its relative motion with respect to the original RF line is determined by maximizing the cross-correlation function. The main deficiency of this method is its inability to account for local compression prior to correlation. A significant compression creates dissimilarities be- tween blocks which causes the cross-correlation method to fail. Accounting for compression by temporal stretching of the compressed signal substan- tially enhances the estimation result of the CC method [16]. However, as both time-shift and temporal stretching must be estimated to maximize the cross-correlation, accurate motion estimation is obtained at the expense of significant computational cost. Another method uses an adaptive stretching algorithm to maximize the correlation between the pre- and post-compression signals [17]. While this method yields very good results for the cases under static compression, the signal-to-noise ratio (SNR) will degrade in the presence of dynamic vibra- tions or when the compression is very small, accounting for limited sensitiv- ity. Conventional CC also suffers from this disadvantage. In determining the quality of a compression estimation algorithm, the SNR is evaluated where the signal is computed as the average value of strain in a region that is supposed to have constant strain, and noise is the standard deviation from that expected value. Usually, the SNR is plotted versus the strain in the medium as a strain filter [18]. For a better statistical assessment, each SNR point in the strain filter can be an averaged SNR value from repeated experiments or different observations. Detecting the zero-crossings in the RF signals has been proposed as another technique to estimate the motion [19]. In that method, the zero- crossings of the original and compressed signals are found and assigned to each other to estimate local tissue motion. Since matching of the zero- crossings is performed by a simple one-to-one assignment procedure, there can be mismatches at areas of high strain, with the result that the strain filter shows a steep drop. It was also shown in [19] that using a multi- compression scheme and averaging the results increases the signal-to-noise ratio of the estimation. A new motion estimation technique called peak search algorithm (PSA) is proposed in this chapter. The main steps of the technique are summarized below: 1. Detect the peaks of the uncompressed and compressed RF signals using the continuous wavelet transform. 18 Chapter 2. Peak Search Strain Imaging 2. For each peak in the compressed signal, find the best match in the uncompressed signal. 3. Define local tissue motion as the matched peak displacement. These steps are explained in detail in Section 2.2. Figure 2.1 shows an RF signal before and after the compression is ap- plied. Arrows indicate the displacements of the peaks which show the tissue motion. It can be seen that due to compression, the peaks preserve their shape while changing their position. The motion of the peaks is caused by local strain in the tissue. So, the motivation of the present method is to track the displacements of the peaks and correlate them with the ax- ial strain in the tissue. This method differs from the zero-crossing method in [19] in two major ways. First, using the wavelet transform to find the peaks of the signals enables the PSA to avoid the peaks that are caused by noise or limited sampling rate. In other words, in this new method not all of the local maxima of the RF signal are considered as signal peaks; rather, by using the wavelet transform, only those that have suitable frequency content would be marked as real peaks. Second, the one-to-one assignment of the zero-crossings is replaced by a more sophisticated procedure to match the peaks in the compressed and uncompressed signals. In summary, a new method is proposed to find and match peaks in pre- and post-compression RF A-lines. The theory of this technique is presented in Section 2.2. Results are shown in Section 2.3, where the algorithm was tested on simulated signals and phantoms. A discussion and conclusion fol- low. Appendix A includes a brief review of wavelet theory and its application for peak detection. A discussion of regularity of the wavelet transform is presented in Appendix B. 2.2 Theory Due to the inhomogeneous nature of real tissue, strain varies along an A-line. Hence, different compression ratios should be considered when estimating the relative motion of two signals. The main assumption made is that tissue can be considered to be locally homogeneous, i.e., due to an applied stress, RF signals preserve their shapes locally while undergoing displacement and compression. To what extent this is true for high strains will be discussed in more detail in Section III. Finding the peaks in the signals robustly is the key aspect of the proce- dure. The peak detection algorithm should not miss the intrinsic RF signal 19 Chapter 2. Peak Search Strain Imaging Time / Distance Am pl itu de Pre−compression RF signal Post−compression RF signal Figure 2.1: A typical RF A-line (solid line) and the same line after com- pression is applied (dotted line). Peak displacements are shown by arrows which give an estimate of the motion inside tissue. The signals are samples taken from an A-line in a homogeneous gelatin phantom, manually com- pressed from the right and imaged using a transducer at the left. No other boundary conditions are imposed. 20 Chapter 2. Peak Search Strain Imaging peaks, while ignoring those caused by high frequency noise or sampling. This aim can not be achieved by the conventional time domain peak search algorithms without pre-filtering of the signals. The accuracy of the peak search procedure is affected significantly by choice of the filter. After signal filtering, either the peaks or the zero-crossings can be sought. Generally, zero-crossings are easier to detect; but peaks can also be found by looking for the zero-crossings in the derivative of the original signal. The advantage of finding the zero-crossings of the derivative is that having differentiated the signal, the regularity of the singular points can also be investigated in order to obtain a more accurate estimate, as is further explained in appen- dix B. The process of filtering and differentiating a signal can be efficiently implemented using the wavelet transform. In similar problems, especially in biomedical signal processing, wavelet methods to find the R-wave in the QRS-complex of the electrocardiogram (ECG) have been successfully used to overcome the interference of high frequency noise and low frequency ar- tifacts [20,21]. The main characteristics of the continuous wavelet transform and its ability to find the peaks and singularities of signals are described in the appendices. Here, it is sufficient to point out that the peaks in the RF signals are detected by finding the zero-crossings in the wavelet domain with a particular choice of the wavelet function. As shown in [22], zero-crossings in the wavelet domain correspond to the peaks in the original signal if the mother wavelet is the first derivative of a Gaussian smoothing function. 2.2.1 Peak Detection Procedure Let Wsf(x) be the wavelet transform (WT) of the RF signal f(x) at time x and scale s, with the mother wavelet ψ(x) derived from a smoothing function θ(x), as: Wsf(x) = f(x) ∗ ψs(x) = s d dx (f ∗ θs) (x) , (2.1) where, ψ(x) = dθ(x) dx , (2.2) θ(x) is a Gaussian distribution with standard deviation σ, θ(x) = 1 σ √ 2π exp [−x2/(2σ2)] , (2.3) 21 Chapter 2. Peak Search Strain Imaging and ∗ represents a convolution. Equation 2.1 shows that the extrema of the smoothed signal f ∗ θs(x) correspond to the zero-crossings in Wsf(x) while the inflection points in the smoothed signal are the extrema of Wsf(x). The terms modulus maxima and minima refer to the maxima and the minima of the continuous wavelet transform. The following steps are in- volved in finding the peaks in an ultrasound RF signal: 1. Calculate the discrete wavelet transform of the signal at scales s1 and s0 (s1 > s0). 2. Find the zero-crossings in the wavelet domain at scale s1, i.e. all points zk for which Ws1f(zk) = 0. 3. Find the modulus maxima and minima of the wavelet transform that correspond to zk in both scales. Using these extrema, calculate the regularity at point zk. 4. If the regularity is positive, proceed; otherwise remove zk from the list of peaks and start over. 5. Reduce the effect of low sampling rate by interpolating the exact lo- cation of the peak (sub-sample approximation). For a 2.5MHz or 5MHz ultrasound transducer, s1 = 4 and s0 = 2 are suitable scale choices, since the pass-band of the filters at these scales have the most overlap with the spectrum of the RF-signals. Appendix A in- cludes a discussion on how to select the proper scales for other values of the ultrasound center frequency. Checking the regularity in step 4 guaranties that only those peaks of the signal are detected that contain the proper fre- quency content. A negative regularity for a stationary point indicates that the signal has more power in scale s0 than s1 in that location. As discussed in Appendix B, this can be a result of high frequency noise interference. Hence, only the peaks with positive regularity values are of interest. Figure 2.2 shows an arbitrary function with its wavelet transform using the first derivative of a Gaussian shaped smoothing function as the mother wavelet. The correspondence between the characteristic points in the signals in Figs. 2(a) and (b) can be observed. However we focus mainly on the peak locations of the original signal which correspond to the zero-crossings of Wsf(x). For a local maximum, the corresponding zero-crossing in the wavelet domain is a transition from negative to positive sign, while for a local minimum it is vice versa. The smoothing function that is used in this 22 Chapter 2. Peak Search Strain Imaging 0 (a) f(x) (b) WT Time/ Distance (c) SAI Time / Distance Time / Distance 0 0 Figure 2.2: Correspondence between the extrema and the inflection points of a signal and zero-crossings and modulus maxima of the wavelet transform: (a) a sample waveform; (b) wavelet transform of the waveform with a mother wavelet that is the first derivative of a Gaussian shaped smoothing function; (c) the mother wavelet used in this example. 23 Chapter 2. Peak Search Strain Imaging study is a Gaussian shaped low-pass filter. The resulting wavelet function is shown in Fig. 2(c). In order to have a more accurate estimation of the peak locations, a sub- sampling technique has been used to interpolate the peak position within a fraction of a sample. By using the value of the signal at the maximum and its two adjacent samples, the peak location can be estimated by interpolation. We have chosen quadratic interpolation for its simplicity, although other interpolation techniques (e.g. cosine) have been used for sub-sample peak estimation in RF signals. After this procedure, the peaks of the RF signal are detected. A vector zk is generated that contains the locations of the peaks. This routine is applied to both pre- and post-compression signals and the peaks are placed in two different vectors. The vectors containing the positions of the peaks in the uncompressed and compressed signals are named zu and zc, respectively. 2.2.2 Peak Matching After the peak vectors have been formed, each peak in zc should be matched to one in zu. Since the deformation is not smooth throughout the tissue, it is not accurate to match the elements in zc and zu one to one. A more sophisticated peak matching algorithm is needed and is described below. For each element in zu, starting from the first, its match in zc is estimated using the average compression ratio of the tissue that was calculated from the previous peak matches. Using the abbreviations CR and EP for the average compression ratio and expected peak location, respectively, we have that EP = CR × zuk . To begin with, the compression ratio is set to 1. Next, zc is searched for the closest peak in the proximity of EP. The element in zc that is the closest to EP and no farther than a pre-specified neighborhood (NB) is selected by finding k such that: min |zck − EP| , subject to, |zck − EP| ≤ NB , NB is specified in number of samples and is originally set to one sample. This means that each peak should be found in the expected location or at most at NB samples away from the expected location. If the search were successful, CR is updated and the next peak is processed. Otherwise, the next peak in zu is examined while NB is increased by one sample. The block diagram for 24 Chapter 2. Peak Search Strain Imaging k = 1 CR = 1 NB = 1 k = k+1, u kzCREP ´= Find m that minimizes EPz cm - subject to, NBEPz c m <- Did such m exist? u kz are c mz matched peaks. NB = NB+1 u k c m z z CRCR aa +-= )1( Figure 2.3: Block diagram of the peak matching procedure. 25 Chapter 2. Peak Search Strain Imaging the peak matching procedure is shown in Fig. 2.3. The compression ratio used at every iteration is the average obtained from previous iterations. α in Fig. 2.3 is a weight that determines the effect of the new estimation on updating the compression ratio. It provides memory to the compression ratio CR and reduces the effect of one false match on the rest of the estimation. In our study α was equal to 0.3. A low value of α ensures that if a false match has occurred at some point in an RF-line, the error will propagate less. The reason is that the average compression ratio (CR) has a higher weight in determining the expected location of the new peak in the compressed signal. On the other hand, a higher value for α makes the algorithm more sensitive to sharp axial strain changes. Any other method, or α value that gives more weight to the previous compression ratio than the new estimation would be suitable for this purpose as well. 2.2.3 Strain Estimation Since the peaks are not linearly distributed in an RF line, the motion along the line can be interpolated from the displacements at the peak locations. Once interpolated, the displacements may then be under-sampled to the desired resolution. The choice of the interpolation method is seen not to have a large effect on the signal-to-noise ratio of the estimation. Four dif- ferent methods of interpolation have been tested including cubic-spline, lin- ear, nearest neighbor and displacement averaging inside overlapping blocks. The differences between the final under-sampled displacements were seen to be 0.1% of the total displacement power. Therefore, instead of interpo- lation and under-sampling, the displacements were simply averaged inside the overlapping blocks for reasons of efficiency. So, once the displacements are estimated, each line is divided into overlapping rectangular windows of specific length and the motion of the peaks that lie inside each block are averaged, giving a single value to the displacement of that block. The re- sulting displacement profile is considered as the displacement along that RF line. To obtain strain, the gradient of the displacement is computed. 2.3 Simulations The algorithm was first tested on a one-dimensional simulated homogeneous tissue under several compression states. The center frequency of the ultra- sound beam was fixed at 2.5MHz and the signals were sampled at 40Mz. The depth of the tissue was 4cm and the imaging window covered 100 par- allel RF lines laterally. RF lines were generated by convolving a Gaussian 26 Chapter 2. Peak Search Strain Imaging shaped point-spread function (PSF) having a center frequency of 2.5MHz with the predetermined scatterer positions. For the compressed signals, the new locations of the scatterers were computed with the known compression factor and the convolution repeated. In order to calculate the estimation signal-to-noise ratio (SNRe) and to analyze the strain filter [23], the homo- geneity is modeled by spreading approximately 20 scatterers per pulse width of the ultrasound beam. Compressions were applied starting from 0.2% up to 20% in steps of 0.2%. This range of strains encompasses the range used in practical strain imaging. The SNRe for each RF line is calculated from the strain estimate in that line. The value of the real strain over the standard deviation in the strain estimate yields the SNRe. For each compression state, the results from 100 lines are averaged to yield a statistical estimate of the accuracy of the algorithm. In order for the peak search algorithm to be suitable for real-time implementation, no prior temporal stretching has been applied to the signals. The accuracy of PSA is compared to a modified cross-correlation method that uses the prior estimates to increase the efficacy and the dynamic range [24]. Based on the motion of the last block, the displacement of the new block is approximated and the correlation is performed in a small region. Sub-sample approximation has been used to estimate the peak of the cross- correlation function more accurately. Temporal stretching would likely in- crease the dynamic range of both methods. Neither the PSA nor the CC was implemented with temporal stretching to improve the SNRe since it would make a large impact on computational efficiency. 2.4 Results 2.4.1 Numerical Results The average value of the estimated strain versus real applied strain is plotted in Fig. 2.4(a). A window length of having a rectangular shape and overlap of 50% is used for this purpose for both methods. Also shown in Fig. 2.4(a) is the variance of the PSA estimator at strains up to 10%. The illustrated curves indicate that PSA is an unbiased estimator for as high as 7% of strain, exhibiting a small variance, especially for lower strains. As can be seen, the dynamic range of the CC method is less than that of the PSA. Strain filters for the same block length and overlap are plotted in Fig. 2.4(b) for different compression states. Under the effective dynamic range of each of these two methods, the SNRe of the PSA is almost three times as high as that of the 27 Chapter 2. Peak Search Strain Imaging 0 2 4 6 8 10 12 14 16 18 200 2 4 6 8 10 12 14 16 18 20 RealStrain (%) Es tim at ed St ra in (% ) Actual strain CC estimator mean PSA estimator variance PSA estimator mean (a) 2 4 6 8 10 12 14 16 18 200 1 2 3 4 5 6 7 8 Strain (%) SN R e PSA CC (b) Figure 2.4: (a) Estimated versus real strain under different compressions for the PSA (circle markers) and the CC method (square markers), superim- posed on the real strain line. The variance of PSA is shown as well for up to 10% of strain (cross markers); (b) Strain filter for the average values of PSA (circles) and CC method (squares). Dotted lines are standard deviations. A window length of 1mm and overlap of 50% has been used for both methods. 28 Chapter 2. Peak Search Strain Imaging CC method as seen in Figure 2.4(b). The large standard deviation of PSA at high frequencies indicates that the peak matching procedure begins to fail and the resulting strain images are not reliable. A high variance is also seen at lower frequencies for the PSA which results from having low density of the scatterers in the RF lines. Note that a higher scatterer density would result in a higher contrast for the peaks of the RF signals, thus likely making PSA more accurate. The effect of window length and overlap on the SNRe is shown in Fig. 2.5. As shown in Fig. 2.5(a), larger segments result in more averaging of the displacement data and obviously increase the SNRe at the expense of lowering the resolution. The overlap of the windows in this case is 50%. The same behavior results from CC and other time domain techniques. Also, for large window lengths, Fig. 2.5(a) shows larger changes of the SNRe vs. strain. This can be due to the fact that by increasing the window length, there are fewer windows in total and hence the averaging may result in less reliable and more fluctuating values. As illustrated in Fig. 2.5(b), lower segment overlaps result in higher values of SNRe. For 1mm window length, three overlap values have been tested. With less overlap, the number of blocks in an RF line decreases and, as a consequence, the strain deviation among blocks drops, and so the SNRe increases. In an estimation problem, the error can be defined in two ways. The first way of measuring the accuracy involves determining the estimation bias in each individual realization. In the strain estimation problem, a small error of this type corresponds to having unbiased estimation in every line of the image. The second way of defining the estimation error is to assume that the estimator is unbiased, then calculating the deviation from the expected value, and finally averaging this value resulted from several observations. Both the strain filter and the real-versus-estimated-strain errors that have been discussed so far and illustrated in Fig. 2.4 are of this second type of error. These errors are only valid in the interval in which the estimator is unbiased. For example, if the average of the estimated strain is not equal to the real strain for a certain compression value, the value of the strain filter at that point is not acceptable. In fact, calculating the SNRe for the strain filter using the average value of the strains as the signal value does not resolve this problem. Likewise, for the error plotted in Fig. 2.4(a), the curves do not yield any information on the strains for the individual lines and how their distributions vary around their expected values. In order to overcome this lack of information on the distribution of the strains of individual lines, we propose a type of root-mean-square (rms) error. If the real value of strain in an experiment is s and there are N lines in which 29 Chapter 2. Peak Search Strain Imaging 2 4 6 8 10 12 14 16 18 200 10 20 30 40 50 60 70 80 Strain (%) SN R e 1mm window 3mm window 5mm window (a) 2 4 6 8 10 12 14 16 18 200 1 2 3 4 5 6 7 8 Strain (%) SN R e 75% overlap 50% overlap 25% overlap (b) Figure 2.5: (a) Strain filters of PSA for different values of (a) window length for 50% of overlap, and (b) segment overlaps having 1mm segment length. 30 Chapter 2. Peak Search Strain Imaging strain has been estimated and its average in line i, (1 ≤ i ≤ N) is ŝi, the rms error can be defined as: Erms = √∑N i=0 (ŝi − s)2 . (2.4) With this definition of error, if there are either many lines in which the average values of the estimated strain have a deviation from the real value of compression or only a few lines that have a significant estimation inaccuracy, the error in equation (2.4) will increase. Fig. 2.6 shows Erms for the PSA and CC method. For both methods the error increases linearly with the real strain. Up to about 7% of strain, the Erms for CC is almost 170 times as high as that of the PSA. A significantly lower error for the PSA compared to the CC method means that the number of lines in which PSA fails in giving a proper estimate is much smaller than that of the CC method. The abrupt increase and fluctuation for the strains higher than 8%, explains that both methods fail to yield an unbiased estimation from that point forward. However, the behavior of the CC method beyond that region depends on the implementation, e.g. on the method of choosing a suitable search region. For the sensitivity analysis, a wider range of compressions has been an- alyzed. Fig. 2.7 shows the SNRe for compressions from 0.001% up to 15%. It can be seen that PSA is not sensitive to strains of less than 0.1%. How- ever even at that point, the SNRe for the PSA is higher than the maximum value of SNRe for the CC method. The reason that PSA is more sensitive to low strains is that CC performs the sub-sample approximation on the cross- correlation function between a few wavelengths of the signals while PSA does the same approximation on the RF peak locations. Although SNRe decreases for high and low strains, these two drop-offs have different causes and hence can be explained differently. For high strains, as the rms error increases, the incorrect estimation will result in a low value for the SNRe. But for low strains, the rms error is low and therefore the estimation has the least bias; however because the signal value is very low, even a small deviation in the estimation reduces the SNRe. Fig. 2.8 shows the effect of different levels of the sonographic noise on the PSA displacement estimation in a 1D simulation. At low strains, the algorithm is sensitive to small additive noises, while at higher strains, the performance is degraded less. Given that an additive noise introduces inaccuracy in estimating the exact location of the peaks, this error, being less than half a wavelength, will be statistically similar for different values of noise. Thus, as expected, at higher compressions the signal to noise ratio is 31 Chapter 2. Peak Search Strain Imaging 2 4 6 8 10 12 14 16 18 2010 -2 10 -1 100 101 10 2 10 3 RealStrain (%) E r m s CC PSA Figure 2.6: RMS errors for the PSA (solid line) and CC method (dashed line). Segment lengths are 1mm and overlap is 50%. 32 Chapter 2. Peak Search Strain Imaging 10-3 10-2 10-1 100 101 102 0 1 2 3 4 5 6 7 Strain(%) SN R e PSA CC Figure 2.7: Strain filters for lower strains to analyze the sensitivity, for PSA and CC method. The solid lines are the averages of the SNRe and the dotted lines depict its standard deviations. 33 Chapter 2. Peak Search Strain Imaging less affected by sonographic noise than at low compressions. For comparison, a 2D phantom with a hard circular inclusion has been simulated and tested with the PSA and CC method. The results are shown in Fig. 2.9. Three states of compression have been tested with each of the methods. It is seen that PSA gives sharper and less noisy elastograms. In this example, PSA starts to show incorrect estimations in a few lines for 7% strain in the bottom-left image of Fig. 2.9. The reason for the errors is the failure of the peak matching procedure. One can see how the error of a few false matches propagates along a line. As strain increases, the peak matching procedure fails, due to the change in the peak shapes. This might be a result of not having enough scatterers, thus having a drastic decrease in the amplitude of the RF signal. Although simple filtering, such as a 2D median filter, could likely remove such artifacts, comparisons are better made on original data prior to filtering. Quantitative comparison of the images of the PSA and CC method is further possible by using the elastographic contrast-to-noise ratio (CNRe) which is defined as follows: CNRe = 2(m1 −m2)2 σ21 + σ 2 2 , (2.5) where m1 and m2 are average strains in the target and the background and σ1 and σ2 are the standard deviations of the strain in those regions [25]. A higher CNRe will result if the estimated strain values in the target and the background have small variances or the contrasts of the average strain in the two regions are low. The values for the images in Fig. 2.9 are 23.9dB, 24.2dB and 10.4dB for the PSA and 9.5dB, 9.2dB and 8.9dB for the CC method for compressions of 1%, 4% and 7% respectively. The corrupted lines in Fig. 2.9(e) under 7% of compression cause the CNRe to drop for PSA. From the simulation results one can conclude that the elastograms resulted from PSA have higher contrast and less noise than those with the CC method but occasional mismatches are possible. Also as discussed in Section 2.5, PSA is computationally more efficient than correlation based methods. 2.4.2 Experimental Results Elastograms for a real 3-layered 3D phantom experiment with the middle layer twice as stiff as the background are shown in Fig. 2.10. 1.6% strain has been applied to the phantom to produce the required deformation. The phantom that was constructed for this experiment was a 3×4×5 cm PVC 34 Chapter 2. Peak Search Strain Imaging 10 20 30 40 50 Inf0 1 2 3 4 5 6 7 Sonographic SNR (dB) SN R e 0.4 % 0.8 % 1.6 % 3.2 % 6.4 % Figure 2.8: Signal-to-noise ratio of the PSA versus sonographic noise level. Different curves represent different levels of compression. Higher strains show smaller performance reduction at high the presence of high sonographic noise, compared to lower strains. 35 Chapter 2. Peak Search Strain Imaging Width (mm) D e pt h (m m ) PSA − 1% Compression 10 20 30 40 10 20 30 0 0.5 1 1.5 2 2.5 (a) Width (mm) D e pt h (m m ) CC − 1% Compression 10 20 30 40 10 20 30 0 0.5 1 1.5 2 2.5 (b) Width (mm) D e pt h (m m ) PSA − 4% Compression 10 20 30 40 10 20 30 2 4 6 8 (c) Width (mm) D e pt h (m m ) CC − 4% Compression 10 20 30 40 10 20 30 2 4 6 8 (d) Width (mm) D e pt h (m m ) PSA − 7% Compression 10 20 30 40 10 20 30 0 10 20 (e) Width (mm) D e pt h (m m ) CC − 7% Compression 10 20 30 40 10 20 30 0 10 20 (f) Figure 2.9: Elastograms of a simulated phantom with a hard inclusion using 1mm widows with 50% overlap. The Top, middle and bottom rows show elastograms for 1%, 4% and 7% compressions using the: PSA (left column), CC method (left column). The scale bars show the range of strain in each image in percent. 36 Chapter 2. Peak Search Strain Imaging phantom held stationary on the bottom and compressed from the top by the imaging transducer. Ultrasound data were captured using the Ultrasonix RP500 (Ultrasonix Medical Corporation, Burnaby, BC) with a 40MHz sam- pling rate and a 5MHz centre frequency. The wavelet transform scales that were used for this experiment are the same as the ones used in the simula- tions. As seen in Fig. 2.10, the harder layer is clearly visible, while it was not possible to delineate it in the B-mode image. Fig. 2.10(a) shows the strain image obtained using PSA and Fig. 2.10(b) shows the strain image resulting from the CC method. The average of the strain in all the lines and the standard deviations for the two methods are shown in Fig. 2.10(c) and 2.10(d). Window length and overlap were 1mm and 50%, respectively, for both methods. It is common practice to filter the noise in the displace- ment estimation, using median filters. Hence, a 5×5 median filter has been applied to the estimated displacements prior to differentiation. While both images are clear elastograms, the CNRe values for the PSA and CC strain images are 11.9dB and 9.3dB respectively. To investigate cases with a larger amount of lateral motion, another experiment was performed using a phantom with a circular inclusion. The stiffness of the inclusion was approximately twice as high as the background. The parameters used to capture RF data were the same ones as used in the previous experiment and the phantom had the same overall dimensions as the three-layered phantom. Fig. 2.11 illustrates the strain images resulting from PSA under four different compression levels. As shown, a very small strain of 0.12% results in lower contrast in the image and poorer visibility of the inclusion compared to larger compressions. This is expected, given the results of Fig. 2.7 where a lower value of SNRe was calculated for smaller compressions. The image quality in Fig. 2.11 is seen to improve when higher compressions up to 1% are applied. A high strain (brighter) region can be seen on the top and bottom of the inclusion in Fig. 2.11(b), 2.11(c) and 2.11(d). This can be simply explained by the fact that the same amount of compression has been applied to all of the lines of the image. Thus, having a low-strain region in any line requires the other regions of that line to have higher than the average strain. 2.5 Discussion A new peak search algorithm was proposed that is aimed to be a low-noise real-time strain imaging technique. The method has been compared with a conventional time domain method that uses the cross correlation of signal 37 Chapter 2. Peak Search Strain Imaging Line Number D e pt h (m m ) PSA Strain Image 20 40 60 10 20 30 1 2 3 (a) Line Number D e pt h (m m ) CC Strain Image 20 40 60 10 20 30 1 2 3 (b) 10 20 300.5 1 1.5 2 2.5 3 3.5 Depth (mm) S tr ai n (% ) PSA strain variation (c) 10 20 300 1 2 3 4 Depth (mm) S tr ai n (% ) CC strain variation (d) Figure 2.10: Elastograms for a real phantom with three layers using 1mm windows and 50% overlap. The middle layer is approximately two times as stiff as the other two layers. (a) Shows the result of the PSA and (b) illustrates the CC result. The Elastograms are computed under 1.6% com- pression. (c) and (d) show the average values of all the lines and the standard deviation for PSA and CC respectively. 38 Chapter 2. Peak Search Strain Imaging Line Number D ep th (m m ) 0.12% Compression 20 40 60 20 40 0.1 0.2 0.3 0.4 0.5 (a) Line Number D ep th (m m ) 0.28% Compression 20 40 60 20 40 0.2 0.4 0.6 0.8 (b) Line Number D ep th (m m ) 0.45% Compression 20 40 60 20 40 0.2 0.4 0.6 0.8 1 1.2 (c) Line Number D ep th (m m ) 1% Compression 20 40 60 20 40 0.5 1 1.5 2 2.5 (d) Figure 2.11: Elastograms for a real phantom with a circular inclusion using 1mm windows and 50% overlap. The inclusion is approximately two times as stiff as the background. Strain images are shown here for cases of having (a) 0.12%, (b) 0.28%, (c) 0.45%, and (d) 1% compression. 39 Chapter 2. Peak Search Strain Imaging segments to estimate the displacement. Compared to previous work, the SNRe values for the CC method in this study seem to be lower, which may be due to several reasons. These include using the gradient instead of a least squares method for the strain calculation, not using temporal stretching, having a lower sampling rate compared to the other implementations and a lower center frequency for the RF signal. The least squares estimation for the strain increases SNRe at the cost of resolution. Temporal stretching will also increase the signal to noise ratio [16], but then the method will be much more difficult to implement in real-time. The same reasoning applies to PSA. As an example, the conventional way of implementing the global stretching involves multiple resampling of the compressed signal at different rates and comparing the correlation results of the resampled signals with the uncompressed RF signal. The resampling and cross-correlating of the entire RF-lines are exhaustive procedures that drastically affect the efficiency and might prevent real-time implementation of the methods. The effect of the sampling rate on SNRe is also discussed in [19]. It has been shown that a higher sampling rate increases the SNRe value. In the present work, compatible with the hardware we use, the sampling frequency has been set to 40MHz. The poor quality of the strain images resulting from the phantom ex- periments are mainly due to the fact that the raw data is being displayed. However, if the images were post-processed using a two-dimensional median filter, the histogram of the images were equalized or simply a least-squares method were utilized to obtain the strain from the displacement field, the SNRe and CNRe values would increase. The improvement in the results for the proposed method compared to the standard cross-correlation technique is due to a different solution strat- egy. The cross-correlation method attempts to find a region in the com- pressed signal that has the maximum similarity to the selected region in the pre-compression signal. However, after compression, the blocks that have maximum correlation with each other may not necessarily correspond to the same location in the tissue, as the local compressions of the RF signals are not taken into account. Therefore using the adaptive stretching as ex- plained in [17] yields great improvements in the results. But in the peak search algorithm, the idea is that after compression the positions of the indi- vidual peaks change but their patterns do not undergo a major deformation. This can be verified by investigating the ultrasound RF lines under different strains and is illustrated in Fig. 2.1. The advantage of using the wavelet transform for finding the peaks is that it filters out the undesired portions of the spectrum and reduces the effect of noise on the peak detection. In fact, 40 Chapter 2. Peak Search Strain Imaging instead of using the wavelet transform, one might be able to use a band-pass filter together with an accurate peak detector to get the same results. How- ever by using the wavelet transform with the mother wavelet explained in Section 2.2.1, these two tasks can be merged into a single convolution with the wavelet function which increases the computational efficiency. In this study, ultrasound center frequencies of 2.5MHz and 5MHz and sampling rate of 40MHz were used. For a different center frequency, the wavelet transform would be computed using the same mother wavelet, but the signals may be analyzed at different scales. Referring to Appendix A, choosing a different scale may make the algorithm more sensitive to high frequency noise. According to Fig. A.1(a), a higher scale would result in de- tection of high frequency peaks in the signals which in turn may increase the sensitivity to noise. One can see in Fig. A.1(a) that for center frequencies of 2.5MHz and 5MHz, a scale of 22 overlaps the best with the signal spectrum. The same discussion on how to choose the scale is valid if the sampling frequency is different from 40MHz. Details on how to choose the proper scales are explained in Appendix A. The effect of the center frequency on the zero-crossing method is similar and is discussed in [19]. The SNRe and CNRe values are supposed to increase for higher ultrasound center frequen- cies, given that the sampling rate increases proportionately. The reason for choosing low center frequencies for the simulations and experiments in the present work is to allow for sufficiently many samples in a wavelength. Given that the ultrasound machine that we use has a fixed sampling fre- quency, using a high center frequency would result in too few samples in a wavelength. This in turn deteriorates the sub-sampling localization of the peaks and reduces the accuracy. However, by adjusting the sampling rate, higher center frequencies can yield higher qualities in the elastograms. The 1D simulation results in Fig. 2.9 show that the SNRe values of PSA are in good coherence with the strain filters presented in Fig. 2.4(b). Compared to the simulations, the experimental results in Fig. 2.10 have lower signal-to-noise ratios. The main reason for this difference is the out- of-plane motion that could not be controlled in the experiment. Also lateral motions have not been compensated for, which results in inaccuracy of the axial motion estimation. An attempt was made to minimize these effects in the experiment by stabilizing the transducer with a robotic arm while compressing the phantom. It should be noted that most of the axial motion estimation techniques, including PSA, are affected by non-axial motion. If the axial and lateral motions are estimated separately, one of the existing lateral tracking methods could be used. The accuracy of PSA is more likely to be affected by the sonographic 41 Chapter 2. Peak Search Strain Imaging noise than that of the correlation based techniques. In the presence of considerable sonographic noise, the location and shapes of the peaks will be corrupted, resulting in possible failure of the peak-finding algorithm. The issue of the sonographic noise in imaging the soft tissue is not as significant as for blood vessels [26]. According to Fig. 2.8, a high signal-to-noise ratio of 30dB, as is the case for imaging soft tissue, makes the effect of sonographic noise small in the peak search method, if the relative strains of the successive frames are higher than 2%. Using wavelet filtering, the robustness of the algorithm to noise (sonographic or other additive noise) is also improved. As shown in Fig. A.1(b), ultrasound RF signals are relatively narrowband signals and any high frequency interference can be filtered to increase the signal-to-noise ratio of the RF signals. Steps 3 and 4 of the peak matching algorithm described in Section 2.2.1 are meant to check the regularity of the wavelet transform at the peak loca- tions. These are optional if the compression is small. For high compressions, the regularity check increases the accuracy. For low strains these steps do not affect the results, and the computational efficiency can be increased by excluding them. In order to have an estimate of the computational efficiency of PSA rel- ative to conventional correlation techniques, one can calculate the number of operations required for both methods. Assuming that there are k blocks, each having N samples, for the case of a 2.5MHz center frequency and 40MHz sampling frequency, under 1% compression, the CC method needs a search region of roughly 50 samples. Therefore, the number of operations required for CC is of the order of [3(2N +50) log2(2N +50)]×k×Cn. Here, Cn is the cost for normalizing the cross-correlation functions and is equal to 3. The first term in the brackets is the cost of convolving the two windows in the compressed and uncompressed signals in the desired search region. The second term in the brackets is due to finding the peak of the cross- correlation function. For the case of the peak search algorithm, the cost would be roughly proportional to 3MW + γP , where M is the number of samples in one line of the RF signal, W is the length of the mother wavelet (equal to 32 for the experiments in this paper) and P is the number of pos- itive peaks in one signal and 3 < γ < 8, depending on the implementation of the algorithm. This term accounts for the cost of finding the peaks in the two signals. The factors that affect γ include: the implementation of the sub-sampling procedure, the choice of sub-sampling approximation in the wavelet-domain or in the time-domain, and the steps taken in order to minimize the impact of false matches or lost peaks. However, the second term is small compared to the cost of the convolution. With the specifica- 42 Chapter 2. Peak Search Strain Imaging tions that have been mentioned for the imaging device, in the case of 1mm segments with 50% overlap, PSA is approximately 12 times faster than the conventional CC method. This speed is sufficient for the PSA to be used in real-time implementations. The method has been written in C++ and run on an Ultrasonix RP500 ultrasound machine with a Pentium IV CPU and 1024MB RAM. A frame-rate of 12fps was achieved for 128-line images, thus approximately 0.7 msec/RF-line. Simplifying the method can increase the frame-rate. As an example, instead of sub-sample approximations of the peak locations, the zero-crossings of the wavelet transform can be accu- rately measured with the sub-sample procedure. Thus, the method would only analyze the wavelet transform of the RF signals which may provide an improvement in the speed of the algorithm. 2.6 Conclusion A novel technique for strain imaging was introduced in this paper that is based on finding the peaks of the ultrasound RF signals. The motivation is that, due to compression, the local signal profile remains approximately the same, but the locations of the peaks change according to the local elasticity of the medium. In order to locate the peaks of the RF signals, a robust peak detection method based on the wavelet transform has been used that can reject portions of the signal outside the spectrum of the RF signal. A new definition for error in strain imaging has also been proposed that evaluates the bias of the estimation in the individual lines of the image. This type of error can be used as complementary information to other accuracy evaluation tools such as the strain filter. Compared to the conventional cross correlation method, the proposed algorithm has a better signal-to-noise ratio, level of error, and speed of processing. It also has the potential for increasing the resolution to the level of inter-peak distance by reducing the window length while maintaining the overlap value. If strain is sufficiently small and the peaks are found and matched correctly, the resolution can be reduced to the level of the ultrasound wavelength for the proposed method, while such a resolution is not achievable with most other methods. Nevertheless, this increase in resolution will only be achieved at the cost of lowering the signal-to-noise ratio. This technique is also a good choice for real-time strain imaging where strain is low and both signal to noise ratio and computational efficiency are important. 43 References [1] M. Bro-Nielsen, “Finite element modeling in surgery simulation,” Pro- ceedings of the IEEE, vol. 86, no. 3, pp. 490–503, March 1998. [2] Y. Zhu, T. J. Hall, and J. Jiang, “A finite-element approach for young’s modulus reconstruction,” IEEE Transactions on Medical Imag- ing, vol. 22, no. 7, pp. 890–901, Sep. 2003. [3] F. Kallel and M. Bertrand, “Tissue elasticity reconstruction using linear perturbation method,” IEEE Transactions on Medical Imaging, vol. 15, no. 3, pp. 299–313, Jun 1996. [4] C. Sumi, A. Suzuki, and K. Nakayama, “Estimation of shear modulus distribution in soft tissue from strain distribution,” IEEE Transactions on Biomedical Engineering, vol. 42, no. 2, pp. 193–202, Feb. 1995. [5] A. R. Skovoroda, S. Y. Emelianov, and M. O’Donnell, “Tissue elasticity reconstruction based on ultrasonic displacement and strain images,” IEEE transactions on ultrasonics, ferroelectrics, and frequency control, vol. 42, no. 4, pp. 747–765, 1995. [6] A. J. Romano, J. J. Shirron, and J. A. Bucaro, “On the noninvasive determination of material parameters from a knowledge of elastic dis- placements: Theory and numerical simulation,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 45, no. 3, pp. 751–759, May 1998. [7] R. Sinkus, J. Lorenzen, D. Schrader, M. Lorenzen, M. Dargatz, and D. Holz, “High-resolution tensor mr elastography for breast tumour detection,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1649– 1664, June 2000. [8] M. Fatemi, A. Manduca, and J. F. Greenleaf, “Imaging elastic prop- erties of biological tissues by low-frequency harmonic vibration,” Pro- ceedings of the IEEE, vol. 91, no. 10, pp. 1503–19, Oct. 2003. 44 Chapter 2. Peak Search Strain Imaging [9] P.-C. Li and M.-J. Chen, “Strain compounding: a new approach for speckle reduction,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 49, no. 1, pp. 39–46, Jan. 2002. [10] R. L. Maurice, J. Ohayon, G. Finet, and G. Cloutier, “Adapting the lagrangian speckle model estimator for endovascular elastography: the- ory and validation with simulated radio-frequency data,” Journal of the Acoustical Society of America, vol. 116, no. 2, pp. 1276–86, Aug. 2004. [11] C. Sumi, “Fine elasticity imaging utilizing the iterative rf-echo phase matching method,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 46, no. 1, pp. 158–66, Jan. 1999. [12] A. Pesavento, A. Lorenz, and H. Ermert, “Phase root seeking and the cramer-rao-lower bound for strain estimation,” in Proceedings of IEEE Ultrasonics Symposium, vol. 2, 1999, pp. 1669–72. [13] A. Pesavento, C. Perrey, M. Krueger, and H. Ermert, “A time-efficient and accurate strain estimation concept for ultrasonic elastography using iterative phase zero estimation,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 46, no. 5, pp. 1057–67, Sep. 1999. [14] M. A. Lubinski, S. Y. Emelianov, and M. O’Donnell, “Speckle tracking methods for ultrasonic elasticity imaging using short-time correlation,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Con- trol, vol. 46, no. 1, pp. 82–96, Jan. 1999. [15] S. G. Foster, P. M. Embree, and W. D. J. O’Brien, “Flow velocity pro- file via time-domain correlation: Error analysis and computer simula- tion,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 37, no. 3, pp. 164–175, May 1990. [16] T. Varghese and J. Ophir, “Enhancement of echo-signal correlation in elastography using temporal stretching.” IEEE Transactions on Ultra- sonics, Ferroelectrics and Frequency Control, vol. 44, no. 1, pp. 173–180, Jan. 1997. [17] S. K. Alam, J. Ophir, and E. E. Konofagou, “An adaptive strain estima- tor for elastography,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 45, no. 2, pp. 461–72, March 1998. 45 Chapter 2. Peak Search Strain Imaging [18] T. Varghese and J. Ophir, “A theoretical framework for performance characterization of elastography: the strain filter.” IEEE Trans Ultra- son Ferroelectr Freq Control, vol. 44, no. 1, pp. 164–172, Jan. 1997. [19] S. Srinivasan and J. Ophir, “A zero-crossing strain estimator for elas- tography,” Ultrasound in Medicine and Biology, vol. 29, no. 2, pp. 227– 238, Feb 2003. [20] C. Li, C. Zheng, and C. Tai, “Detection of ecg characteristic points us- ing wavelet transforms,” IEEE Transactions on Biomedical Engineer- ing, vol. 42, no. 1, pp. 21–28, Jan. 1995. [21] H. Eskandari and M. B. Shamsollahi, “Detection of the characteristic points in electrocardiogram using the time-scale transforms,” in Pro- ceedings of the International Congress in Biological and Medical Engi- neering (ICBME), Dec. 2002. [22] S. Mallat and S. Zhong, “Characterization of signals from multiscale edges,” IEEE Transactions on Pattern Analysis and Machine Intelli- gence, vol. 14, no. 7, pp. 710–32, Jul 1992. [23] T. Varghese, J. Ophir, E. Konofagou, F. Kallel, and R. Righetti, “Tradeoffs in elastographic imaging,” Ultrasonic Imaging, vol. 23, no. 4, pp. 216–248, Oct. 2001. [24] R. Zahiri-Azar and S. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, Oct. 2006. [25] T. Varghese and J. Ophir, “An analysis of elastographic contrast-to- noise ratio.” Ultrasound in Medicine and Biology, vol. 24, no. 6, pp. 915–924, July 1998. [26] E. Konofagou and J. Ophir, “A new elastographic method for esti- mation and imaging of lateral displacements, lateral strains, corrected axial strains and poisson’s ratios in tissues,” Ultrasound in Medicine and Biology, vol. 24, no. 8, pp. 1183–1199, Oct 1998. [27] S. Mallat and W. L. Hwang, “Singularity detection and processing with wavelets,” IEEE Transactions on Information Theory, vol. 38, no. 2, pp. 617–643, Mar 1992. 46 Chapter 3 Viscoelastic Parameter Estimation Based on Spectral Analysis 2 3.1 Introduction Mechanical properties of soft tissue are known to be affected by malignancies and abnormalities. Stiffness is the mechanical property that has been most widely studied. Previous papers have shown that benign breast lesions are four times as stiff as the background tissue and carcinomas are eight times stiffer than fibroadenoma [1]. Ultrasound elastography has emerged as a new imaging modality for depicting variations in tissue stiffness, but sen- sitivity and specificity are still being explored. For example, in a recent paper, fourteen patients with breast pathology were studied, six of whom appeared to have benign lesions after biopsy and seven had malignancies [2]. Ultrasound elastography, however, could not discriminate the benign tissue, but was successful in detecting most of the malignant tumors. Imaging of the mechanical properties in soft tissue to investigate or monitor the patho- logical changes in the body, generally referred to as elastography, has been widely studied in the past two decades [1–12]. While elasticity correlates with the physical structure of soft tissue, vis- cosity can be a means to describe metabolism-dependent features. Carcino- mas with high blood vessel concentration exhibit higher viscosity than the surrounding area. In [13], fifteen patients with different breast pathologies were examined and the cancerous breast tissue was seen to be on average 4.4±3.2 times more viscous than the surrounding tissue. Estimation of the tissue viscosity has recently been studied by a number of groups. Chen et al. analyzed the dispersion of narrow-band shear waves in homogeneous soft tis- 2A version of this chapter has been published. H. Eskandari, S. E. Salcudean, and R. Rohling, “Viscoelastic Parameter Estimation Based on Spectral Analysis”. IEEE Trans. on Ultras. Ferr. Freq. Cont., vol. 55, no. 7, pp. 1611-1625, Jul. 2008 47 Chapter 3. Viscoelastic Parameter Estimation sue [14]. From the phase of the dissipating wave at different locations, they could estimate the shear wave velocity which, in turn, is related to the shear elasticity and viscosity. Catheline et al. measured the overall viscoelastic properties by analyzing the speed and attenuation of a transverse wave at several frequencies [15]. Insana et al. used ultrasound to image the deforma- tion of soft tissue due to a compressional step force [16]. Using the theory of viscoelastic solids, they modeled the relaxation of the material by two exponentials: one with a low and the other with a high time-constant. At frame-rates less than 4 frames/second, they were able to detect the dynamic creep and relaxation in the tissue which enabled them to quantify a small time-constant of around 1s and a long time-constant of nearly 200s [17]. Bercoff et al. measured the temporal variation of the shear wave inside a homogeneous phantom and compared it with the waveform of a simulated viscoelastic model [18]. The parameters that would yield the most simi- larity between the simulation and the experimental waveforms were related to the mechanical parameters of the homogeneous phantom. Sinkus et al. described a technique to reconstruct the viscosity as well as the anisotropic components of elasticity using MR elastography [19]. While they obtained a good separation between malignant, benign and background breast tissue using elasticity, the shear viscosity did not seem to be a significant and dis- criminating feature. Acoustic remote palpation introduced in [9] has been utilized by Girnyk et al. to create a dynamic strain field and to estimate the dynamic properties of the tissue [20]. They were able to characterize the phantoms and tissue samples in vitro at the focal point of vibration, assuming homogeneity in the medium. In order to measure the dynamic properties of the vitreous body in the human eye, Walker et al. quantified the step response of the tissue due to a remotely induced acoustic radiation step force [21]. Work is ongoing to improve sensitivity and specificity of these techniques to make them suitable for routine clinical use. As an alternative approach, Turgay et al. showed a new technique that modeled soft tissue with a series of Voigt elements and added inertia [11]. The proposed reconstruction technique involved solving an inverse problem for the model parameters such that the error between modeled and mea- sured displacements was minimized. Preliminary simulation results showed promise for the reconstruction of the viscosity as well as elasticity. Also, spectral analysis of the time series of displacements was utilized to yield a robust elasticity estimation of the medium. The transfer functions calcu- lated between displacements at different locations were used to characterize the viscoelastic properties independent of the model assumption. Salcud- ean et al. used the concept of transfer functions to delineate the prostate 48 Chapter 3. Viscoelastic Parameter Estimation boundaries in phantom and experiments in vivo [12]. The transfer function concept will be explored in more detail in the present work. In particular, the phase and the imaginary parts of the trans- fer functions will be shown to be suitable means for identifying the dynamic properties of the tissue. The relationships between viscoelastic properties and transfer functions will be explored through simulations and phantom experiments and validated using a rheometer and small samples of material. The overall purpose is to study the fundamentals of viscoelastic parameter identification in a controlled environment. In order to characterize the ability of imaging methods to measure vis- coelastic properties of tissue there is a need for tissue mimicking materials with controllable static and dynamic properties. There is a broad collection of literature on elasticity phantoms using gelatin mixtures [22,23], agar com- posites [24] or polyvinyl alcohol (PVA) [25]. Various techniques have been employed to induce viscosity contrast in phantoms. It is known that increas- ing the concentration of gelatin in the gelatin-based phantoms will increase both the elasticity and the viscosity of the material [26]. At low concentra- tions (i.e. less than 10%) the dependency of the elasticity parameter on the gelatin concentration was seen to be higher than that of the viscosity. There- fore, increasing the proportion of gelatin in water can be used to mimic a hard and viscous lesion with a higher time-constant [17,21]. However, based on the rheology data shown in section 3.5, at concentrations above 10%, the time-constant of the gelatin phantom remains nearly unchanged with respect to the gelatin percentage in water. Therefore, increasing the gela- tine concentration results in a stiffer and more viscous phantom while the relaxation-time of the material does not change. Using agents such as glyc- erol [11, 20], glutaraldehyde [27] or formaldehyde [17, 23] that increase the amount of cross-linking in the gelatin is another way to elevate the viscous effect. In most cases, the cross-linking will affect the elastic and viscous properties as well as the speed of sound in the tissue. In this work, we in- duce a local inclusion in the phantom by using a material different from the background. It is shown that PVA sponge can be used as an inclusion in a gelatine based phantom to control the contrast of the static or dynamic mechanical parameters within the tissue. 49 Chapter 3. Viscoelastic Parameter Estimation Figure 3.1: One dimensional model of soft tissue with springs, dampers and masses, assembled as Voigt elements. In the present work, the tissue is considered bound at one end and excited dynamically from the other end. By measuring the displacements with the help of ultrasound, the parameters of the model can be identified. 50 Chapter 3. Viscoelastic Parameter Estimation As an overview, this paper aims to reconstruct the local viscoelastic parameters of soft tissue. The reconstruction techniques are verified by sim- ulating a one-dimensional (1D) network of Voigt elements as a mechanical model of tissue. A transient force is generated and applied to the tissue model and the displacements at designated locations are measured. Using the acquired displacements, the tissue parameters are identified as explained in the following sections. The ability of the proposed methods to detect elas- tic and viscous lesions in soft tissue are tested on tissue mimicking materials. The experimental setup is capable of vibrating a small sample with a filtered wide-band excitation, measuring the force, and observing the motion ultra- sonically over a short period of time. Next, the magnitude and phase of the displacement spectra are analyzed and associated with the viscoelastic properties of the material using well-known parameter identification tech- niques. From the resulting property calculations, cross-sectional images of viscoelastic properties are generated. 3.2 Models Soft tissue is known to have nonlinear behavior due to external excitations. Effects such as hysteresis, large deformation and nonlinear and time-varying mechanical properties can generate higher order harmonics. These nonlin- earities are seen to have rather insignificant effect on the response of most soft tissues at small strains [27]; thus, it is viable to adopt a linear framework in order to analyze and model the response of the soft tissue. 3.2.1 Linear Viscoelastic Model A mechanical structure comprises elastic properties, viscous properties and inherent density. The elastic behavior of a structure makes it resilient to an external force or deformation. A purely elastic element will vibrate in- finitely in response to a momentary excitation. The elasticity in soft tissue can be modeled by a finite number of spring elements in the direction of compression, assuming the compression to be unidirectional. However, due to internal friction and material relaxation effects, a damping term is intro- duced in the vibration. This effect, which is a liquid-like behavior of soft tissue, is mostly induced by the viscosity of the material. In the theory of viscoelastic solids, viscosity can be modeled by damper elements. Assum- ing linearity, the force-velocity relationship of the dampers is linear, as well as the force-displacement relationship of the springs. Lastly, any element 51 Chapter 3. Viscoelastic Parameter Estimation that bears mass exhibits a linear relationship between force and accelera- tion. Hence, tissue can be modeled as a set of interconnected elements, each being composed of a spring, a damper and a point mass. Usually the inertia in the soft tissue can be ignored at low frequencies, considering the relatively low density of human tissue, which is close to that of water. It will be shown in Section 3.5 that this approximation is valid in the experiments at frequencies below 50Hz, depending on the phantom or tissue material. Therefore at low frequencies, the tissue can be modeled by a series of springs and dampers. For an incremental region in the tissue, these two components can either be in parallel, which results in the Voigt model, or in series, which gives the Maxwell model. A one-dimensional model consisting of several Voigt elements with added inertia is shown in Fig. 3.1. It has been shown that the Voigt model is more accurate in describing the dynamics of the soft tissue [15]. The constitutive relation between the stress (σ) and strain (ǫ) of a Voigt element is a linear first order differential equation as follows: Eǫ(t) + ηǫ̇(t) = σ(t) , (3.1) where E is the Young’s modulus or elasticity and η is the viscosity. By ignoring the densities in Fig. 3.1, the model will contain a series of Voigt elements [28]. The equation of motion for this model can be obtained by discretizing eqn (3.1). For a typical element i, this gives: ki[ui(t)− ui+1(t)] + bi[u̇i(t)− u̇i+1(t)] = f(t) , (3.2) where ki is the stiffness of the i th element and bi is the viscous damping coefficient. f(t) is the applied force and is the same at every point in the network, since the effect of inertia has been ignored. The interested reader can refer to [11] for a detailed analysis of such a model as well as the rela- tionships between the discrete model parameters and the tissue mechanical properties. When a dynamic compression is applied to a non-homogeneous sam- ple of finite dimensions, the resulting motion will be determined by several factors. These factors include the method of excitation, the boundary condi- tions, mode conversion, reflection, refraction, nonlinear response and other complex phenomena. This investigation attempts to create axial motion and minimize other motions through careful design of the experimental appara- tus. As will be shown in Section V, by applying a distributed force on the surface of the phantom, the out-of-plane components of the stress tensor can 52 Chapter 3. Viscoelastic Parameter Estimation be minimized and a state of plane-stress can be approximated. Boundary conditions and excitation method are also carefully controlled. The result is that the axial motion and its spatial variations can be tracked by means of ultrasound and a tissue model can be used to identify the mechanical parameters in the tissue. Two key assumptions are made for modeling the tissue response: (1) the axial strain response of tissue to the axial excitation is linear; (2) the linear relationship between axial strains at different tissue locations in response to the axial excitation can be described by a first order lead-lag relationship, i.e., an elastic component and a viscous component arranged as a Kelvin- Voigt model [11, 21, 28–30]. Such a dynamic system can account for the energy dissipation in the system and justify the observed phase change. The first assumption can be supported by the measurement of the coherence function between the excitation and the estimated displacements. This value which is seen to be higher than 0.96 in the experiments in Section 3.5, provides a measure of the proportion of the input energy that is transferred into axial tissue motion in a linear fashion [31]. Introducing a comprehensive model that predicts all the nonlinear ef- fects present in the medium is beyond the scope of this paper. In fact, in the aforementioned linear 1D model, the lateral motions and non-symmetric deformations are ignored. To justify the use of such a model in experiments, a uniform and symmetric excitation should be applied to the phantom. Also, the boundary conditions should be controlled to minimize the shear effects and induce a low-speed compressional wave in the medium. When these con- ditions are satisfied and a nearly plane-stress state can be maintained in the medium, the speed of the longitudinal wave from the excitation (in contrast to approximately 1540m/s for ultrasound), can be as low as c = √ E/ρ for finite cross-sectional volumes with free sides [32]. In soft tissues, the Young’s modulus (E) is approximately 10kPa and the density (ρ) is close to 1000kg/m3, which result in a phase velocity of c = 3.2m/s. These calcula- tions are valid for wavelengths that are large compared to the cross-sectional dimension of the sample [33]. Since pulse-based ultrasound is used to mea- sure the tissue motion, it does not provide an instantaneous snapshot of the displacements at a specific time, but given the relatively small excitation frequency (<30Hz), slow excitation phase velocity (3.2m/s) and relatively fast pulse velocity (1540m/s), this source of error is negligible. Assuming linearity, a given excitation frequency will produce localized tissue motion at that frequency at all locations within the tissue, so the pulse repetition frequency of ultrasound is set greater than twice the maximum applied ex- citation frequency to satisfy the Nyquist requirement. As will be shown in 53 Chapter 3. Viscoelastic Parameter Estimation Section V, by applying proper boundary conditions and excitation, the axial motion of the tissue from the excitation can be easily tracked by means of ultrasound and the proposed model can be used to identify the mechanical parameters in the tissue. The experiments are constructed in such a way that the axial motion is the most dominant component in the medium so the model is appropriate. Achieving or circumventing these same conditions in experiments in vivo is an acknowledged challenge and is the subject of ongoing study. 3.2.2 Spectral Analysis and the Transfer Functions If the linear network in Fig. 3.1 is excited by a single frequency force wave- form, the displacements and local strains will have the same frequency, but different amplitudes and phase lags. If the force spectrum is wide-band, the spectra of the displacements and strains will also cover a range of frequen- cies. The shape of each spectrum depends on the model parameters, i.e. the mechanical properties of tissue. The transfer functions between the applied force and the nodal displace- ments were studied in previous work [11]. Wherever the force could not be measured, the transfer functions were computed by taking the displace- ment of one of the blocks as the reference. The low-frequency asymptotes of the transfer functions were correlated with the quasi-static properties of the tissue to delineate the inclusions and the anatomical boundaries [12]. Appendix C includes a review of the transfer function calculation between the stress and the strain/displacement waveforms. Estimation of the elastic modulus from the displacement transfer functions is described in [11] and is repeated in Appendix D using the strain transfer functions. 3.2.3 Power-Law Effect It is known that human soft tissue and most tissue mimicking phantoms behave like non-Newtonian fluids in the presence of a dynamic excitation. Newtonian fluids, like water, preserve their elastic and viscous properties at all frequencies, while non-Newtonian materials exhibit different elastic and viscous properties at different frequencies. In most materials, the variation in viscosity versus frequency can be modeled by a power-law [34]. For a generalized power-law fluid, the viscosity can be expressed as: η = ηc ω n , (3.3) 54 Chapter 3. Viscoelastic Parameter Estimation where ηc is the consistency index, ω is the angular frequency and n is the flow index. For Newtonian fluids, n = 0. The case of n > 0 is rarely encountered in biological tissues and is referred to as dilatency or the shear-thickening effect. The case of n < 0 is known as the shear-thinning or pseudo-plastic effect and is very common in polymers and gelatin composites. Accounting for the power-law effect in viscosity is a special case of the Kelvin-Voigt fractional derivative model of the viscoelasticity, which has been previously used to model the response of canine tissue over a wide frequency range [29]. 3.3 Algorithms The magnitude and phase of the transfer functions depend on the elastic and viscous properties of the medium. To show how transfer functions between the strains and the applied stress change as an inclusion appears in the medium, a four-element network of mass-spring-dampers is simulated. The second element is assumed to have higher elasticity and the third element is assumed to have higher viscosity than the background. The density is assumed to be constant and equal to 1000kg/m3. Figures 3.2(a) and 3.2(b) show the values of elasticity and viscosity in the simulated medium. The model was excited for 10 seconds with 40Hz low-pass filtered white Gaussian noise as the force, and the displacements and strains were calculated. The transfer functions between the applied stress and the strains at every element were then computed and their magnitudes and phases are plotted in Figs. 3.2(c) and 3.2(d). The transfer functions were computed using blocks of 1s with 75% overlap. Each curve in these two figures corresponds to the transfer function for one of the elements in the system. One can observe the change in the low-frequency asymptotes of the magnitudes due to different elastic properties. The phase change, however, follows a more subtle trend which can be associated with the ratio between the viscosity and elasticity. This ratio will be called the tissue relaxation-time and can be reconstructed from the phase or imaginary part of the transfer functions. For the second element with a smaller relaxation-time, the phase is closer to zero, while for the third element with a larger relaxation-time, the phase of the transfer function is larger. A small spatial drift in the phase at low frequencies results from the influence of inertia. Second order effects can also be seen in the transfer functions, causing the peaks in the magnitude plot and large variations in the phase plot at frequencies higher than 20Hz. 55 Chapter 3. Viscoelastic Parameter Estimation 1 2 3 4 10 20 30 40 50 Element Index El as tic ity (K Pa ) Actual Elasticity (a) 1 2 3 4 10 20 30 40 50 Element Index Vi sc os ity (P a.s ) Actual Viscosity (b) 100 101 0 2 4 6 8 x 10−3 Frequency (Hz) |Hε σ (ω )| TF Magnitude Element 1 Element 2 Element 3 Element 4 (c) 100 101 −0.4 −0.3 −0.2 −0.1 0 ∠ H ε σ (ω ) Frequency (Hz) TF Phase Element 1 Element 2 Element 3 Element 4 (d) Figure 3.2: The effect of changing the viscoelastic parameters on the transfer functions. (a) and (b) show the elasticity and the viscosity of a four-element one-dimensional simulated phantom. (c) and (d) are the magnitude and phase of the transfer functions (TF) between the strains at the applied stress. The excitation force was white noise, low-pass filtered at 40Hz. 56 Chapter 3. Viscoelastic Parameter Estimation 3.3.1 Parameter Estimation If the frequency is sufficiently small, the tissue motion will be governed by only the elastic properties and the dynamic effects do not play a significant role. This is also evident from Fig. 3.2(c) where the magnitudes of the transfer functions at low frequencies appear to be nearly flat. Thus, one can obtain an estimate for the elasticity of the medium from the low-frequency asymptotes of the transfer functions [11,12], (see also Appendix D). The phase change in the transfer function is due to both elasticity and viscosity variations. As seen in Fig. 3.2, the phases of the transfer functions change considerably when an inclusion occurs in either the viscosity or the elasticity. It is shown in Appendix C that the transfer function between the local strain and the applied stress can be parameterized in terms of the Voigt model parameters as follows: Hǫσ(ω) = 1 (E + jωη) . (3.4) Assuming ωη is much smaller than E, it follows that: ωη E ≈ −∠Hǫσ(ω) . (3.5) The above equation shows that a high viscosity region affects the phase in the same way as a low elasticity region. While one can estimate the viscosity by calculating η from eqn (3.5), it is more straightforward to define and estimate a relaxation-time (τ) for viscoelastic tissue as follows: τ = η E ≈ − 1 ω ∠Hǫσ(ω) . (3.6) A large τ indicates a slow response for the corresponding element. This is equivalent to the group delay of the transfer function at that frequency, considering that the phase changes linearly versus frequency. If the viscosity has a power-law trend and elasticity is independent of fre- quency, the relaxation-time will be subject to a similar behavior. The distri- bution of the flow index resulting from the power-law in the relaxation-time will thus yield another tissue parameter to be identified. If this power-law trend is ignored, the relaxation-time can be assumed constant in a frequency range that inertia does not interfere. Therefore, instead of dealing with a single frequency, the slope of the phase with respect to frequency in eqn (3.6) can be calculated to obtain a robust estimate of the relaxation-time. If the stress is unknown, the transfer functions can be computed with 57 Chapter 3. Viscoelastic Parameter Estimation respect to the strain at a specific location. Hence, equation (C.5) from Appendix C can be used to derive the following estimate for the distribution of relaxation-time: τ̂ = τ0 − 1 ω ∠Hǫǫ0(ω) , (3.7) where τ0 = η0/E0 is the relaxation-time for the element of reference. Since τ0 is unknown, it can be eliminated from eqn (3.7) by assuming it to be constant. Although this is a biased estimate, the bias is constant for all of the elements and is equal to τ0. Therefore, the above equation gives a relative but correct profile of the distribution of the relaxation-time in the tissue. However, the power-law effect cannot be readily seen here without measuring the force. If the power-law behavior is ignored, τ̂ can be averaged at several frequencies to obtain more robust estimates. In order to quantify the estimation error, the stress is assumed to be known in the simulations in section 3.4. However, the experimental data in section 3.5 do not include any force measurement, thus relative estimates are presented. 3.4 Simulations In order to investigate the accuracy and noise-sensitivity of the proposed methods, one-dimensional finite element simulations were performed. The system in Fig. 3.1 was constructed with 100 elements. A band-limited force was applied and the displacements were computed. The details on how to simulate a mass-spring-damper network are presented in Appendix E. The model was assumed to represent a cross-section of 1mm2 and a depth of 50mm of the tissue. The medium was homogeneous with the elasticity of 20kPa and viscosity of 20Pa·s. These values were chosen arbitrarily in the range of the viscoelastic properties for average soft tissues. The density of the phantom was 1000kg/m3 which is equal to the density of water and most soft tissues. The simulation frame-rate was 1kHz and the duration of each test was 10 seconds. The applied stress was white noise low-pass filtered at 30Hz. The simulated model was linear, hence the magnitude of the applied force would not cause any nonlinear effects in the displacements, nor would it change the estimation results. The maximum strain that was observed in the phantom was no more than 1%. The displacements were computed at each node and the transfer functions were calculated between the strains and the applied stress. Using eqns (D.1) and (3.7), relative estimates for elasticity and relaxation- 58 Chapter 3. Viscoelastic Parameter Estimation time were calculated. For the simulated homogeneous phantom, the para- meter estimates for all of the elements are expected to be equal; thus a high standard deviation indicates poor estimation. If p0 is the actual value of the parameter and the vector p̂ is the estimate for all of the elements, the estimation error (Ere) can be defined from the standard deviation (STD) as follows: Ere = STD(p̂) p0 . (3.8) The sensitivity of the estimations to the noise level can be measured by calculating Ere under the assumption that the displacements are noisy. Therefore, the model was simulated and the actual strains were calculated. A certain level of noise was then added to the strains and the transfer functions were calculated using the noisy strain data. This noise is meant to account for the quantization error and limited resolution in the displacement estimation, interference of the lateral motion of the tissue in the axial motion tracking, non-linearities in the tissue, and imperfections of the model. A white Gaussian noise model has been assumed to address the sum of these sources of error. A more detailed study of the displacement and strain noise distribution is beyond the scope of this paper, but has been investigated by others [35]. A simple noise model has been adopted in order to compare the effects of various parameters on the accuracy of the estimation, rather than precise quantification of the error. The accuracies of the proposed methods are evaluated for a wide range of strain signal-to-noise ratios (SNR). For each level of the added noise on the strain data, one-hundred simulations were performed. The accuracy of the elasticity and relaxation-time estimations were calculated each time using eqn (3.8) and the results were averaged for all of the one-hundred instances. Fig. 3.3 shows the Ere for the elastic modulus and the relaxation-time. The elasticity was reconstructed from the low frequency asymptote of the transfer functions magnitudes. The relaxation-time was calculated from the phase, averaged between 2–10Hz as described in section 3.3.1. The estimation error for the relaxation-time has been divided by ten in order to clearly display both curves in the same figure. The strains in the previous calculations were computed by taking the derivative of the displacements. Alternatively, one can use the least-squares (LSQ) technique to calculate the gradient of the displacement which will reduce estimation variance at the cost of reduced contrast [36]. The effect of the LSQ kernel size on the estimation error is illustrated in Fig. 3.4(a). The network was simulated with 100 elements. The first 50 elements had a back- 59 Chapter 3. Viscoelastic Parameter Estimation 101 10−3 10−2 10−1 100 Strain SNR (dB) Er e Elasticity Estimation Error 10−1 × Relaxation−Time Error Figure 3.3: The error in estimating the elasticity (star markers) and the relaxation-time (RT, circle markers) for different levels of signal-to-noise ratio (SNR) of the strain. The error has been calculated based on the deviation of the estimates from the actual value of that parameter. Note that the estimation error for the relaxation-time is divided by ten. 60 Chapter 3. Viscoelastic Parameter Estimation ground elasticity and relaxation-time equal to 20kPa and 2ms respectively, while the values for the second 50 elements in the network were 40kPa and 1ms. The strain was assumed to have a finite SNR of 15dB. The LSQ kernel size was increased from 2 points, which is the discrete derivative, to 15 points which corresponds to 7.5mm in the tissue. For as many as 100 tests, Ere was calculated for the first half and the second half of the region separately and the results were averaged. It can be seen that increasing the LSQ kernel size reduces the estimation variance for both elasticity and relaxation-time. Note that the increase in the average background parameters compared to the results in Fig. 3.3, marginally enhances the estimation accuracy. Fig. 3.4(b) shows the effect of the LSQ kernel size on the contrast-to-noise ra- tio (CNR) for the two-layered simulated phantom. The CNR is defined as follows: CNR = 2(M1 −M2)2 STD1 2 + STD2 2 , (3.9) where M1 and M2 are the means and STD1 and STD2 are the standard deviations of the estimated parameter at the two regions [35]. A large dif- ference in mean values or very low noise levels can result in high CNR. The CNR resulting from the relaxation-time is multiplied by ten to show both curves in the same figure. In the experimental studies that will be presented in the following section, an LSQ filter with an appropriate kernel size is used to achieve better noise rejection. Note that for these tests, the frequency ranges used for the estimations were kept unchanged to maintain consis- tency. However, in practice, finding a proper frequency range will be crucial to obtaining good parameter estimates and therefore contrast. 3.5 Experiments The proposed methods were tested on tissue mimicking phantoms. Special hardware and software were implemented for controlled tissue excitation and data acquisition. The experimental set-up, phantom construction technique and the experimental results are described. 3.5.1 Actuator and Rheometer Design The schematic of the experimental set-up is shown in Fig. 3.5. The design consists of an immobile platform that seats the actuation system, sensors and the transducers. The motion of the vibrating stage is guided by a linear bearing and preloaded by a spring from the bottom to maintain contact be- 61 Chapter 3. Viscoelastic Parameter Estimation 1 3 5 7 9 11 13 15 10−3 10−2 10−1 LSQ Kernel Size Er e Elasticity Estimation Error 10−1 × Relaxation−Time Error (a) 1 3 5 7 9 11 13 15 101 102 103 104 LSQ Kernel Size Es tim at io n CN R LSQ Effect on CNR Elasticity CNR 10 × Relaxation−Time CNR (b) Figure 3.4: The effect of the least-squares (LSQ) kernel size on the esti- mation. These are calculated using a two-layered simulated phantom. (a) The estimation error Ere for elasticity and relaxation-time decrease as a wider LSQ kernel is applied. (b) The effect of the LSQ kernel size on the contrast-to-noise ratio (CNR). 62 Chapter 3. Viscoelastic Parameter Estimation tween its wheel and a motor-driven tumbler. The linear bearing is actuated by a coreless DC motor (#118778, Maxon Motor, Sachseln, Switzerland) with a ”tumbler” - a cylinder mounted obliquely on its shaft. The angled cylinder acts like a continuously variable elliptical cam. The amplitude is manually adjusted by moving the motor in the transverse direction. An in- cremental encoder with 1000 cycles per revolution (Type E2 from US Digital, Vancouver, WA, USA) was mounted on the motor shaft for position feed- back. The motion of the vibrating stage with respect to the static section is monitored by a position sensing device (PSD) that tracks the light spot of a small infra-red light-emitting diode (IR-LED). The PSD signal was used to validate the motion tracking methods and to control the average strain in the phantom. A one-dimensional PSD (type S3932, Hamamatsu Photon- ics K.K., Hamamatsu City, Japan) together with the required circuitry was mounted on the static part and the IR-LED was connected to the vibrating stage. The control of the motor, as well as the required data acquisition, was done on a PC using a PCI interface card and a terminal board (Q8, Quanser Inc., Markham, ON, Canada). The specifications of the device are summarized in Table 3.1. A real-time program was developed using WinCon (by Quanser) and Simulink (Mathworks Inc., Natick, MA, USA) to generate the desired wave- forms. A velocity-control feedback loop was established between the signal recorded from the encoder and the waveform applied to the motor. The user can either specify a frequency of rotation, which would apply continuous full rotary motion to the motor, or generate any desired wide-band excitation, which would cause the tumbler to oscillate without making complete turns. In the former case, the transverse position of the tumbler determines the vibration amplitude while in the latter case, the power of the generated waveform determines the amplitude of vibration. The device can either be used as a rheometer or as an ultrasound vibro- elastography system. When used as a rheometer, a compression load-cell (Transducer Techniques, Temecula, CA, USA) is mounted on the stage and a small material sample is placed between the vibrating plate and the load- cell. An oscillation with a specified frequency and amplitude is applied and the force and displacement are recorded from the load-cell and the PSD respectively. Precise knowledge of the actual dimensions of the specimen is necessary to calculate the exact values of the elasticity and viscosity. However, the relaxation-time can be determined independently from the sample dimensions and can be recovered from the phase difference between the force and the position signal using eqn (3.7). Note that the stress and the force have the same phase and there is only a scaling factor between 63 Chapter 3. Viscoelastic Parameter Estimation them equal to the cross-section of the sample. Also note that the position recorded from the PSD can be divided by the length of the sample to obtain the average strain. Therefore, the transfer function phase in eqn (3.7) can be simply substituted by the phase of the transfer function between the force and the PSD position signal. When the device is used for vibro-elastography, a linear array ultrasound transducer is mounted in an aluminum case and attached to the stage. A phantom is placed on the vibrating plate and is pressed against the ultra- sound probe on top of it, using ultrasound coupling gel. In this way, a desired excitation can be applied to the motor and the motion within the phantom can be monitored using ultrasound motion estimation techniques. Sequences of radio frequency (RF) data frames can be recorded using a Sonix RP ul- trasound machine (Ultrasonix Medical Corp., Richmond, BC, Canada). For the experiments reported, the recorded ultrasound RF lines have a center frequency of 5MHz and were sampled at 20MHz. After the acquisition, the A-lines were up-sampled by a factor of three to increase the accuracy in measuring the displacements. The acquisition frame-rate depends on the number of lateral A-lines in the imaging window. The axial displacements are estimated using a time-domain cross-correlation technique with prior estimates [37]. The correlation coefficients were monitored to ensure accu- rate motion tracking. The frame of reference for motion tracking was only changed when the average of the correlation coefficients had a significant drop (i.e. below 85%). Thereby, the displacement estimation drift was min- imized and higher accuracy was achieved at low frequencies. Depending on the experiment, the kernel size and overlap can be adjusted; however for a typical experiment that will be described in the next section, they were held constant at 1.5mm and 40%, respectively. 3.5.2 Phantom Construction and Rheometry To construct a phantom with a significant contrast in the relaxation-time, a small piece of polyvinyl alcohol (PVA) sponge (Ceiba Technologies, Chan- dler, AZ, USA) was embedded in a cubic gelatin phantom. This type of sponge exhibits a significantly different viscosity than the gelatin. One can investigate the difference between the relaxation-times in the sponge and gelatin through simple palpation by observing the recovery times. To mea- sure the exact relaxation-times, a piece of sponge and a soft and a hard sam- ple of gelatin phantom were examined using the aforementioned rheometer. The hard gelatin phantom was approximately three times harder than the soft one. Figure 3.6(a) shows the relaxation-times versus frequency. The 64 Chapter 3. Viscoelastic Parameter Estimation Figure 3.5: The schematic of the experimental setup. The ultrasound probe is held stationary, while a DC motor moves a stage and causes the phantom to be compressed. 65 Chapter 3. Viscoelastic Parameter Estimation 3 10 30 10−4 10−3 10−2 10−1 100 Frequency (Hz) R el ax at io n− Ti m e (s) Phantom Rheometry Gelatin−A Gelatin−B Sponge−A Sponge−B Sponge−C (a) 3 10 30 100 101 102 103 Frequency (Hz) Er ro r ( %) Rheometry Error for Gelatin and Sponge Error in Sponge−A Rheometry Error in Gelatin−A Rheometry (b) Figure 3.6: (a) Relaxation-times of the PVA sponge and soft and hard gelatin phantoms versus frequency. Gelatin-A and Gelatin-B are soft and hard gelatin phantoms respectively. Sponge-A is the squeezed PVA sponge after it has been soaked in water, Sponge-B is the partially squeezed sponge and Sponge-C is the wet PVA sponge instantly after being removed from the water. The relaxation-times for the soft and hard gelatin are almost the same while that of the sponge is higher by about one order of magnitude. (b) shows the errors in measuring the relaxation-times of the sponge and gelatin with respect to the data from a stress-controlled rheometer. 66 Chapter 3. Viscoelastic Parameter Estimation sample dimensions were the same, having 3mm thickness and 20mm × 12mm surface area. The relaxation-time of the PVA sponge was measured in three different states. First, the test was performed instantly after the sponge was removed from the water. Then, the sponge was partially compressed and the test was repeated. Finally, complete compression was applied to the sponge to remove the existing water drops and another measurement was taken. The force and the position were recorded and the relaxation-time of the material was estimated as described in section 3.5.1. It can be seen in Fig. 3.6(a) that the two gelatin samples have approximately the same relaxation-times, while that of the sponge is significantly higher. The decreasing trend in the curves with respect to the frequency is due to the power-law effect. These results have also been validated by shear rheometry using a commer- cial stress-controlled rheometer (Bohlin C-VOR, Malvern Instruments Ltd., Worcestershire, UK). Using the generalized Hook’s law and the assumption of incompressibility of the material, it can be shown that rheology with both rheometers should yield the same value for the relaxation-time. Figure 3.6(b) shows the relative error of the measured relaxation-times with respect to the results from the Malvern rheometer. Compared to the rheometry validation results, the curves in Fig. 3.6(a) have 20% and 12% deviation for the soft gelatin and wet sponge respectively. Note that the experiments with both rheometers have variability that depends on experimental conditions which may account for the higher error at high frequencies in Fig. 3.6(b). The gelatin phantom was made using bovine skin gelatin (type G9382- B, Sigma-Aldrich Inc., Oakville, ON, Canada) in water and 2% cellulose (by weight) as the scatterers. The phantom dimensions were 65mm axially, 55mm laterally and 35mm elevationally. The sponge sample, having been soaked in water and compressed periodically over the course of several weeks to remove all air bubbles, was embedded in the gelatin immediately after being removed from the water. The speed of sound in the wet sponge was approximately 1600ms−1 which is close to that of water. However, the attenuation of the ultrasound in the sponge was considerably higher than that in gelatin which can be explained by the high viscosity of its constituent material. Hence, a thin sample was cut for the experiments to ensure good quality in the RF A-lines without shadowing. For fully developed ultrasound speckle, the amplitude distribution of each RF frame should have a Rayleigh distribution with a ratio of mean over standard deviation equal to 1.9 [38]. For the case where the phantom was held stationary, this parameter was 1.91 for all of the frames. This indicates that the recorded A-lines had fully developed speckle distributions, so standard tracking algorithms are used. 67 Chapter 3. Viscoelastic Parameter Estimation 0.9 0.92 0.94 0.96 0.98 1 100 102 104 106 108 Correlation Coefficient Nu mb er o f O ccu ren ce Histogram of CC (a) 100 101 102 10−5 100 Frequency (Hz) No rm aliz ed Ma gni tud e Spectrum Nodal Displacement PSD Measurement (b) 100 101 102 0 0.2 0.4 0.6 0.8 1 Frequency (Hz) Co her enc e V alu e Coherence Function (c) 100 101 102 0 1 2 3 4 5 Frequency (Hz) |H(ω )| Transfer Function −0.1 0 0.1 0.2 0.3 0.4 ∠ H (ω) [rad ] Magnitude Phase (d) 0 10 20 30 40 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Depth (mm) Re lati ve Ela stic ity Estimated Elasticity Elasticity (1−Coherence) Inclusion Boundaries (e) 0 10 20 30 40 −6 −4 −2 0 2 Depth (mm) Re lati ve Re lax atio n−T ime (ms ) Estimated Relaxation−Time 5Hz 15Hz 25Hz Inclusion Boundaries (f) Figure 3.7: Single-line data acquisition and parameter estimation in a gelatin phantom with a PVA sponge inclusion and excited at 3–30Hz. (a) Histogram of the correlation coefficients. (b) Normalized spectrum of the displacement of a node at 30mm (black) and the spectrum of the measurement from the PSD (gray). The coherence function between these two signals is shown in (c). The magnitude and phase of a sample transfer function between element 18 and 30 are shown in (d). The reference element 30 is located at a depth of 17mm. (e) shows the elasticity estimates from the asymptotic magnitude of the transfer functions. The values are normalized with respect to the elasticity of the reference element. The error in estimation is depicted as one minus the value of the coherence function at each point. (f) shows the relaxation-time estimates from the phase of the transfer functions at three different frequencies. Again, the values are the difference between the actual relaxation-time of the phantom and that of the reference element at 17mm. The inclusion is delineated in both figures. 68 Chapter 3. Viscoelastic Parameter Estimation 3.5.3 Data Analysis and Results The transfer functions are computed for all of the strains on each line with respect to the strain of an arbitrary block in the middle of that line. The coherence functions are also calculated to ensure the linearity of the tissue response. A high coherence function - i.e. close to one - indicates that the input-output relation is linear and the noise is small [39]. If the applied stress is known, the coherence function between the stress σ(t) and the strain ǫ(t) at some location in the phantom is: Cǫσ(ω) = |Pǫσ(ω)|2 Pσσ(ω)Pǫǫ(ω) . (3.10) In the presence of noise and non-linearities, the coherence values likely de- crease as the input power level increases [40]. In elastography experiments, this may be due to the larger lateral motion in the phantom or sources of nonlinearities in the tissue such as the effect of large strain and hysteresis. In the first test, the gelatin part of the phantom contained 18 wt.% gelatin powder. A high gelatin concentration was used in order to reduce the elasticity contrast between the sponge and the gelatin. A small piece of PVA sponge (5mm thick) was placed in the middle of the gelatin phantom, approximately 13mm away from the probe face. The phantom was excited with Gaussian white-noise, band-limited to 3–30Hz. The standard deviation of the Gaussian distribution was 80µm and the average strain applied to the phantom was no more than 0.3%. RF-lines were collected for 60 seconds at 1300 frames/second in the Doppler mode from a single A-line passing through the sponge. The imaging window length was 40mm extending from the probe face. The displacements were estimated with 70 blocks of 1.5mm each with 40% overlap between them. The accuracy of the displacement estimates can be verified by examining the correlation coefficients [41]. The histogram of the correlation coefficients at all of the blocks for the duration of the experiment is shown in Fig. 3.7(a), exhibiting an average of 99.4%. To ensure faster motion tracking, stretching has not been applied to the compressed RF signals. However, the high correlation coefficients obtained indicate that the accuracy was well preserved in the results. The spectrum of the displacement of a typical block at 30mm from the probe is shown in Fig. 3.7(b). The spectrum of the PSD measurement is also shown on the same plot. The coherence function between this nodal displacement and the PSD signal is shown in Fig. 3.7(c). The high coherence within the excitation frequency range indicates that the estimated displacements accurately fol- low the position measured by the PSD. The strain was calculated at every 69 Chapter 3. Viscoelastic Parameter Estimation time-step with a 7-point LSQ filter which corresponds to 4mm in the tissue. The transfer and coherence functions for the strains were calculated using blocks of 1 second with 50% of overlap. The reference block to compute the transfer functions was element 30, in the middle of the phantom, with high correlation coefficients in its vicinity. Any other node with high correlation coefficients is an acceptable choice as the reference. The magnitude and the phase of the transfer function for element 18 at 10mm depth are shown in Fig. 3.7(d). Note that the transfer functions are only acceptable at 3–30Hz. Beyond this frequency range, the coherence functions drop and the spectral estimations would be inconsistent. The coherence functions between the lo- cal strains and the strain at the reference block had an average of 0.98 at excitation frequency range. This indicates the high linearity of the tissue response. Higher coherence values result for the nodes closer to the probe, due to the drop in ultrasound power and resolution with depth and lower motion tracking accuracy for deeper nodes. The estimation of the elastic- ity, using the asymptotic magnitude of the transfer functions, is shown in Fig. 3.7(e) and the relaxation-time estimates from the phase of the transfer functions at three different frequencies are plotted in Fig. 3.7(f). It can be seen in Fig. 3.6(a) and 3.7(f) that the contrast between the relaxation-times of gelatin and sponge decreases as the frequency is increased. The error in estimation is shown as one minus the value of the coherence function at each point. A low coherence function indicates non-linearity of the response or high amount of noise in the data, while a coherence value of nearly one, indicates the validity of the linear model. Note that the elasticity estimates are normalized with respect to the unknown elasticity of the reference ele- ment (element 30 at a depth of 17mm) and the estimated relaxation-time of each element is the actual relaxation-time minus the relaxation-time of the reference element, which is also unknown. In another experiment, a soft gelatin phantom (12 wt.% gelatin in wa- ter) was constructed with two inclusions inside it. One inclusion was a hard cylinder of the same gelatin material (18 wt.% gelatin in water), approx- imately 2.5 times stiffer than the background. The other inclusion was a small piece of PVA sponge, saturated with water. Multiple A-lines were captured over the region of interest at 98 frames/second for a duration of 20 seconds. The dimensions of the phantom, excitation characteristics and the displacement estimation parameters were the same as in the previous experiment. The imaging window was 40mm×40mm with 64 RF-lines in the lateral direction. The 1D displacements were estimated with 70 blocks for all the lines, resulting in an average of 98.6% for the correlation coefficients. The strains were calculated with a 4mm LSQ filter. The transfer functions 70 Chapter 3. Viscoelastic Parameter Estimation Width (mm) D ep th (m m) B−mode Image 0 10 20 30 40 0 10 20 30 40 (a) Width (mm) D ep th (m m) Relative Elasticity 0 10 20 30 40 0 10 20 30 1 1.5 2 2.5 3 (b) D ep th (m m) Width (mm) Relative Relaxation−Time (ms) 0 10 20 30 40 0 10 20 30 −1 0 1 2 3 (c) Figure 3.8: Estimated parameters for a gelatin phantom (12 wt.% in water) with a hard gelatin (18 wt.% in water) inclusion at the middle left and a PVA sponge inclusion at the middle right. The inclusions are delineated in the B-mode image (a). The inclusions are not easily detectable in the B-mode image while they are known to be within the dashed ellipses. The elasticity estimate of the phantom can be used to detect both inclusions (b). The low-frequency asymptotes of the transfer functions between 3Hz and 10Hz were used. Note that both inclusions appear as harder regions in the elas- ticity image. While elastic modulus does not provide a contrast between the hard gelatin and sponge inclusions, the relaxation-time estimate in (c) can distinguish the sponge in the gelatin environment. The phase of the transfer functions at 15Hz was used to calculate relaxation-time. As expected from rheometry, the relaxation-times of the soft and hard gelatin materials are nearly the same while sponge has a significantly different relaxation-time. 71 Chapter 3. Viscoelastic Parameter Estimation were calculated for the strain of each line with respect to an element close to the probe (element 5), using the same parameters as the previous exper- iment. The coherence values had an average of 0.9 at the frequency range of excitation. In order to remove possible outliers and present a clear image of the estimated parameters, a 5×5 pixels median filter has been applied to the images. The low-frequency asymptotes of the transfer functions between 3–10Hz were used in order to estimate the relative value of the elasticity in the phantom. Fig. 3.8(b) shows the estimated elasticity, which clearly de- lineates a hard inclusion at the middle left and another one at the middle right of the phantom. The left-side inclusion is the hard gelatin while the right-side inclusion is the sponge. The relaxation-time was estimated from the phase of the transfer functions at 15Hz and is shown in Fig. 3.8(c). The difference between the relaxation-times of the gelatin and the sponge is apparent in this image. 3.6 Discussion The simulations and phantom experiments demonstrate the feasibility of re- constructing the relaxation-time in order to distinguish different soft tissues. The elasticity and the relaxation-time were both estimated using a frequency analysis of the tissue response. Although using a transfer function approach eliminates the need for a specific tissue model, a linear viscoelastic model has been used to verify the algorithms. The rheometry experiments reveal the difference between the relaxation- times for the gelatin and the PVA sponge. From Fig. 3.7(f), the difference between the relaxation-times of the gelatin and sponge is smaller at 15Hz and 25Hz than at 5Hz in Fig. 3.7(f). The decrease in the relaxation-time is mainly governed by a power-law which is evident in Fig. 3.6(a). However, unless the force is measured or the parameter τ0 (the relaxation-time of the reference element) in eqn (3.7) is known, this power-law cannot be deci- phered with the current approach. Table 3.2 contains the consistency index (ηc/E) and the flow index (n) of the best power-law fit for each relaxation- time curve extracted from the data in Fig. 3.6(a). The consistency index explains the higher relaxation-time of the sponge compared to the gelatin and the index n describes the higher negative slope of the curves pertain- ing to the PVA sponge. Both parameters, if successfully estimated, could be suitable features for distinguishing different materials in phantoms and likely in human tissue. Note that here, the elasticity is assumed to be independent of frequency. In this way, the flow index of the viscosity would be the same 72 Chapter 3. Viscoelastic Parameter Estimation as that of the relaxation-time. However, it is shown in [34] for some animal tissue that the elasticity increases slightly as the frequency is increased. As a result, the flow index for elasticity, viscosity and relaxation-time would be different. According to the rheometry results in Fig. 3.6(a), the difference between the relaxation-times of the sponge and the gelatin at 5Hz, 15Hz and 25Hz are approximately 11ms, 4.3ms and 2.4ms respectively. This difference was estimated to be 7.7ms, 3ms and 1.8ms respectively at those frequencies based on the estimation results depicted in Fig. 3.7(f). Considering that using an LSQ filter in the estimations reduces the contrast, the parameter estima- tion and the rheometry data seem to match reasonably well. Figure 3.8(c) estimates this difference in relaxation-time to be approximately 3ms, which also agrees with the rheometry. The presence of tissue harmonics at around 30Hz in the simulations had a large effect on Ere. While the system was simulated by second order mass-spring-damper elements, identification was performed by fitting a first order Voigt model to each element. As a result, the error should be larger for higher frequencies of excitation. This problem was not encountered in the experiments. A first order system could be fit accurately to the transfer functions and the effect of inertia could be ignored at the frequencies below 30Hz. In the experiments, a band-limited force has been applied to the phan- tom. Fig. 3.7(b) shows that most of the excitation power is in the range of 3–30Hz. However, a significant 60Hz interference can be seen in the recorded signal from the PSD which is due to the cross-talk with the power lines and the oscillations in the ambient light. The limited resolution of the PSD accounts for the noise level at higher frequencies. There is also a notice- able peak in the spectrum of the estimated displacement at 86Hz and a few smaller peaks at higher frequencies. This peak, which was not in the exci- tation and is not present in the PSD spectrum, is likely from the harmonic vibrations of the tissue. This identifies the frequency at which the inertia starts to interfere in the spectrum. Some excitation power is also seen in the range of 30–60Hz which is due to the particular shape of the tumbler, producing a nonlinear transform between the rotation of the motor shaft to the motion of the stage. As a result, Fig. 3.7(c) shows coherence values of around 0.7 at this frequency range, but this range was not used for analysis. In the single A-line experiment, the relaxation-time seems to have higher contrast and a smoother profile compared to the elasticity. This is expected as the phantom was designed to have a stiffness close to that of the PVA sponge. The relatively high elastic modulus of the gelatin may have caused 73 Chapter 3. Viscoelastic Parameter Estimation Parameter Value 3dB Frequency Bandwidth 66 Hz Maximum Motion Range 4.42 mm PSD Linearity Range 5 mm PSD Noise Level (std) 9 µm Maximum Measurable Force 2.5 lbf Table 3.1: Rheometer Specifications Material ηc/E [ms] n Gelatin-A 4.6 -0.66 Gelatin-B 4.5 -0.65 Sponge-A 29.2 -0.73 Sponge-B 43.9 -0.84 Sponge-C 62.9 -0.94 Table 3.2: Power-law fit of relaxation-time for the PVA sponge and gelatin 74 Chapter 3. Viscoelastic Parameter Estimation a non-homogeneous sol-gel transition at the time the phantom was made, which can account for the presence of large local variations in the elasticity of the phantom. The other reason for the high variance of the elasticity and the relaxation-time at the depth interval of 25–40mm in Fig. 3.7 is the lower correlation coefficients resulting from the displacement estimation at that section of the phantom. The average correlation coefficients for the first and second halves of the phantom are 99.9% and 98.9% respectively. Considering the high frame-rate of 1,300Hz and the small compression be- tween consecutive frames, the second half of the phantom exhibits relatively small correlation coefficients. Higher lateral motion, poorer quality of the RF A-lines as a result of high attenuation in the sponge, or higher axial motion and signal decorrelation in the phantom might be the sources of the degraded motion tracking at that region. In Fig. 3.4(a), it is shown that using an LSQ filter improves the estima- tion for both elasticity and relaxation-time. This was expected, since the least-squares filter gives a smooth approximation of the derivative at the cost of contrast. It raises the question of what is the appropriate LSQ filter size. The contrast of the elasticity estimate drops by increasing the LSQ kernel size, but the relaxation-time estimation does not seem to suffer from this trade-off. In the range that is shown in Fig. 3.4, the relaxation-time CNR increases as the filter size is increased, so it may be desirable to increase the LSQ filter length to achieve lower estimation error and higher CNR for the relaxation-time. As a trade-off, increasing the filter length limits the resolution and lowers the ability to detect small inclusions. It should also be noted that an extremely wide LSQ filter eventually decreases the CNR of the relaxation-time estimates. While quantitative analysis of a 1D experiment is depicted in Fig. 3.7, the same technique has been extended to a 2D image in Fig. 3.8, using several 1D estimations. Therefore, the parameters have been separately estimated in each line and then the images were formed by joining the indi- vidual lines together. Figure 3.8 shows that the elasticity of both inclusions are higher than the background while the relaxation-time is higher only for the sponge. These results show that the relaxation-time, as a dynamic pa- rameter in the tissue, can be used as a complementary feature to elasticity to distinguish between different soft tissues. In an ideal case, all the components of the displacement tensor and their spatial derivatives are used to reconstruct the three dimensional distribution of the mechanical parameters in tissue. As proposed recently [42], the de- formations due to a shear excitation can be more accurately reconstructed in three dimensions (3D) by imaging the tissue motion ultrasonically from 75 Chapter 3. Viscoelastic Parameter Estimation several angles. While such a method, being capable of solving the wave equation in 3D, can be very promising, it requires the subject to remain still for a few minutes during the examination. The estimation of the parame- ters (especially the dynamic ones) will also require accurate synchronization between the imaging transducer and the excitation at different angles. The simplicity that the linear 1D model in the present work offers makes possible the real-time viscoelastic parameter estimation in soft tissue with a simple excitation scheme. The shortcomings of a 1D model of deformation can be compensated by its advantages in terms of a short period of examination, high computational speed and relative ease of implementation. 3.7 Conclusions It is feasible to estimate the viscoelastic properties of a material by vibrating it over a range of low frequencies and analyzing the time series of motion measurements at multiple locations. This can be performed by calculation and analysis of transfer functions in the frequency domain. From the phase of the transfer functions, the property of relaxation-time can be measured. The relaxation-time is proposed as a potential complementary feature to elasticity to distinguish materials with different mechanical properties. The errors of the algorithm have been quantified through computer simulations. The apparatus for testing these concepts on tissue-mimicking phantoms was described and used to investigate the key factors affecting accuracy. The relaxation-times of different phantom materials were measured successfully. The accuracy of the measured relaxation-times was found to be dependent on appropriate filtering in the strain calculations, suitable frequency range of excitation, and an understanding of the non-Newtonian behavior of the material. The viscosity of the material was modeled with the power-law and the model parameters determined by rheometry. The model parameters calculated from rheometry were similar to the parameters calculated with the transfer functions. Overall, the transfer function method shows promise for future work on measuring the viscoelastic properties of human tissue. 76 References [1] J. Ophir, I. Cespedes, B. Garra, H. Ponnekanti, Y. Huang, and N. Mak- lad, “Elastography: ultrasonic imaging of tissue strain and elastic mod- ulus in vivo,” European Journal of Ultrasound, vol. 3, no. 1, pp. 49–70, Jan. 1996. [2] C. Pellot-Barakat, M. Sridhar, K. K. Lindfors, and M. F. Insana, “Ul- trasonic elasticity imaging as a tool for breast cancer diagnosis and research,” Current Medical Imaging Reviews, vol. 2, no. 1, pp. 157–164, Feb. 2006. [3] R. M. Lerner, K. J. Parker, J. Holen, R. Gramiak, and R. C. Waag, “Sono-elasticity: medical elasticity images derived from ultrasound sig- nals in mechanically vibrated targets,” in Acoustical Imaging. Vol.19. Proceedings of the 16th International Symposium, 1988, pp. 317–27. [4] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastog- raphy: a quantitative method for imaging the elasticity of biological tissues,” Ultrasonic Imaging, vol. 13, no. 2, pp. 111–34, April 1991. [5] L. Gao, K. J. Parker, S. K. Alam, and R. M. Lerner, “Sonoelasticity imaging: theory and experimental verification,” Journal of the Acousti- cal Society of America, vol. 97, no. 6, pp. 3875–86, June 1995. [6] A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov, “Shear wave elasticity imaging: a new ultrasonic tech- nology of medical diagnostics,” Ultrasound in Medicine and Biology, vol. 24, no. 9, pp. 1419–1435, 1998. [7] M. Fatemi and J. F. Greenleaf, “Vibro-acoustography: An imaging modality based on ultrasound-stimulated acoustic emission,” Proceed- ings of the National Academy of Science, U.S.A., vol. 96, no. 12, pp. 6603–6608, June 1999. [8] J. Ophir, S. K. Alam, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, and T. Varghese, “Elastography: ultrasonic estimation and imaging 77 Chapter 3. Viscoelastic Parameter Estimation of the elastic properties of tissues,” Proceedings of the Institution of Mechanical Engineers, Part H (Journal of Engineering in Medicine), vol. 213, no. H3, pp. 203–33, 1999. [9] K. Nightingale, M. Palmeri, R. Nightingale, and G. Trahey, “On the feasibility of remote palpation using acoustic radiation force,” Journal of the Acoustical Society of America, vol. 110, no. 1, pp. 625–634, 2001. [10] J. Bercoff, M. Tanter, and M. Fink, “Supersonic shear imaging: A new technique for soft tissue elasticity mapping,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 51, no. 4, pp. 396–409, 2004. [11] E. Turgay, S. E. Salcudean, and R. Rohling, “Identifying the mechan- ical properties of tissue by ultrasound strain imaging,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 221–235, 2006. [12] S. E. Salcudean, D. French, S. Bachmann, R. Zahiri-Azar, X. Wen, and W. J. Morris, “Viscoelasticity modeling of the prostate region us- ing vibro-elastography.” in Medical Image Computing and Computer- Assisted Intervention - MICCAI 2006. 9th International Conference. Proceedings, Part I (Lecture Notes in Computer Science Vol. 4190). Copenhagen, Denmark: Springer, 2006, pp. 389–96. [13] R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink, “Viscoelastic shear properties of in vivo breast lesions measured by MR elastography,” Magnetic Resonance Imaging, vol. 23, no. 2, pp. 159–65, 2005. [14] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity and viscosity from measurement of shear wave speed dispersion,” Journal of the Acoustical Society of America, vol. 115, no. 6, pp. 2781–2785, 2004. [15] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelka- ram, and J. Culioli, “Measurement of viscoelastic properties of homo- geneous soft solid using transient elastography: An inverse problem ap- proach,” Journal of the Acoustical Society of America, vol. 116, no. 6, pp. 3734–3741, 2004. [16] M. F. Insana, C. Pellot-Barakat, M. Sridhar, and K. K. Lindfors, “Vis- coelastic imaging of breast tumor microenvironment with ultrasound,” 78 Chapter 3. Viscoelastic Parameter Estimation Journal of Mammary Gland Biology and Neoplasia, vol. 9, no. 4, pp. 393–404, Oct. 2004. [17] M. Sridhar, J. Liu, and M. F. Insana, “Viscoelasticity imaging using ultrasound: parameters and error analysis.” Physics in Medicine and Biology, vol. 52, no. 9, pp. 2425–2443, May 2007. [18] J. Bercoff, M. Tanter, M. Muller, and M. Fink, “The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 51, no. 11, pp. 1523–1536, 2004. [19] R. Sinkus, M. Tanter, S. Catheline, J. Lorenzen, C. Kuhl, E. Son- dermann, and M. Fink, “Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography.” Magnetic Reso- nance in Medicine, vol. 53, no. 2, pp. 372–387, Feb 2005. [20] S. Girnyk, A. Barannik, E. Barannik, V. Tovstiak, A. Marusenko, and V. Volokhov, “The estimation of elasticity and viscosity of soft tissues in vitro using the data of remote acoustic palpation,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 211–219, 2006. [21] W. F. Walker, F. J. Fernandez, and L. A. Negron, “A method of imag- ing viscoelastic parameters with acoustic radiation force.” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1437–1447, Jun 2000. [22] T. J. Hall, M. Bilgen, M. F. Insana, and T. A. Krouskop, “Phantom materials for elastography,” IEEE Transactions on Ultrasonics, Ferro- electrics and Frequency Control, vol. 44, no. 6, pp. 1355–1365, 1997. [23] E. L. Madsen, M. A. Hobson, H. Shi, T. Varghese, and G. R. Frank, “Tissue-mimicking agar/gelatin materials for use in heterogeneous elas- tography phantoms,” Physics in Medicine and Biology, vol. 50, no. 23, pp. 5597–5618, 2005. [24] R. Q. Erkamp, S. Y. Emelianov, A. R. Skovoroda, and M. O’Donnell, “Nonlinear elasticity imaging: theory and phantom study.” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 51, no. 5, pp. 532–539, May 2004. [25] K. C. Chu and B. K. Rutt, “Polyvinyl alcohol cryogel: an ideal phan- tom material for MR studies of arterial flow and elasticity.” Magnetic Resonance in Medicine, vol. 37, no. 2, pp. 314–319, Feb 1997. 79 Chapter 3. Viscoelastic Parameter Estimation [26] C. Joly-Duhamel, D. Hellio, and M. Djabourov, “All gelatin networks: 1. Biodiversity and physical chemistry,” Langmuir, vol. 18, no. 19, pp. 7208–7217, 2002. [27] R. Sinkus, J. Bercoff, M. Tanter, J.-L. Gennisson, C. El-Khoury, V. Ser- vois, A. Tardivon, and M. Fink, “Nonlinear viscoelastic properties of tissue assessed by ultrasound.” IEEE Trans Ultrason Ferroelectr Freq Control, vol. 53, no. 11, pp. 2009–2018, Nov 2006. [28] Y. C. Fung, Biomechanics: Mechanical Properties of Living Tissues. Springer, 1993. [29] M. Kiss, T. Varghese, and T. Hall, “Viscoelastic characterization of in vitro canine tissue,” Physics in Medicine and Biology, vol. 49, no. 18, pp. 4207–18, Sept. 2004. [30] T. G. Mezger, The Rheology Handbook, 2nd ed. William Andrew Publishing, 2006. [31] J. S. Bendat and A. G. Piersol, Engineering Applications of Correlation and Spectral Analysis, 2nd ed. Wiley-Interscience, 1993. [32] L. M. Brekhovskikh and V. V. Goncharov, Mechanics of Continua and Wave Dynamics. Springer, July 1985. [33] H. Kolsky, Stress Waves in Solids. Dover Publications Inc., 1963. [34] N. Phan-Thien, S. Nasseri, and L. Bilston, “Oscillatory squeezing flow of a biological material,” Rheologica Acta, vol. 39, no. 4, pp. 409–17, 2000. [35] T. Varghese and J. Ophir, “An analysis of elastographic contrast-to- noise ratio.” Ultrasound in Medicine and Biology, vol. 24, no. 6, pp. 915–924, Jul 1998. [36] F. Kallel and J. Ophir, “A least-squares strain estimator for elastogra- phy.” Ultrasonic Imaging, vol. 19, no. 3, pp. 195–208, Jul 1997. [37] R. Zahiri-Azar and S. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, Oct. 2006. 80 Chapter 3. Viscoelastic Parameter Estimation [38] R. Wagner, M. Insana, and S. Smith, “Fundamental correlation lengths of coherent speckle in medical ultrasonic images,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 35, no. 1, pp. 34–44, Jan. 1988. [39] J. S. Bendat, Nonlinear System Techniques and Applications. John Wiley & Sons Inc., Wiley-Interscience Publication, 1998. [40] V. Z. Marmarelis, “Coherence and apparent transfer function measure- ments for nonlinear physiological systems,” Annals of Biomedical En- gineering, vol. 16, no. 1, pp. 143–157, 1988. [41] S. Srinivasan, F. Kallel, and J. Ophir, “Estimating the elastographic signal-to-noise ratio using correlation coefficients.” Ultrasound in Medi- cine and Biology, vol. 28, no. 3, pp. 359–368, Mar 2002. [42] 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 tomography for in vivo breast examination,” in Proceedings of the IEEE Ultrasonics Symposium, Vancouver, BC, Canada, 2006, pp. 1037–1040. 81 Chapter 4 Viscoelastic Characterization of Soft Tissue from Dynamic Finite Element Models 3 4.1 Introduction Malignant tumors and carcinomas have different mechanical properties com- pared with normal tissue. For many years, physicians have located tumors in soft human tissues by palpating the patient. In the last two decades, many studies have revealed the significance of tissue elasticity (e.g. [1, 2]) and viscosity (e.g. [3,4]) in classifying normal, benign or malignant masses. Several methods have been developed to estimate the local parameters in the viscoelastic medium of soft tissues. Generally, via an imaging modal- ity such as ultrasound or magnetic resonance imaging (MRI), the internal motion of the body is estimated in response to an excitation and can be analyzed to reconstruct the local variations in the viscoelastic properties of tissue. Surveys of such methods can be found in [1, 2, 5, 6, 8]. Ultrasound and MRI are the most popular imaging modalities for track- ing tissue deformation due to an excitation. The reconstruction of the elas- ticity from the boundary force and the internal displacements measured with ultrasound elastography can be formulated in several ways. While a three dimensional (3D) modeling scheme is the most realistic, traditional 2D ultrasound imaging (the most common approach) has limitations on the accuracy of measuring the three components of motion. Axial displacements are measured with the highest accuracy, lateral displacement measurement is less accurate, and elevational displacement measurement is often impossi- ble (note that in ultrasound and ultrasound elastography, axial is along the 3A version of this chapter has been published. H. Eskandari, S. E. Salcudean, R. Rohling and J. Ohayon, “Viscoelastic Characterization of Soft Tissue from Dynamic Finite Element Models”. Physics in Medicine and Biology, vol. 53, no. 22, pp. 6569–6590, Nov. 2008. 82 Chapter 4. Dynamic Finite Element Inverse Problem direction of the excitation and pulse propagation, lateral is perpendicular to axial but within the imaging plane, and elevational is perpendicular to the imaging plane.). These limitations have given rise to techniques that sim- plify the model or ignore one or two components of the displacements. If we consider only the axial displacement and a one-dimensional (1D) model, the reciprocal of the quasi-static axial strain can be interpreted as the local elas- ticity. The viscosity or the relaxation-time can be estimated using a model that predicts the displacement or strain as a function of frequency [6,8–10]. The moduli reconstructed using 1D models suffer from artifacts, because the effect of the boundary conditions on strain is not decoupled from the effect of material properties. To account for axial and lateral displacements, the partial differential equations that describe the stress-strain relationship should be evaluated in a plane. Assumptions on the deformation profile such as being in a plane-stress or plane-strain state can help simplify the general 3D equations. Furthermore, assuming a nearly incompressible and isotropic medium signif- icantly reduces the number of elasticity parameters and leaves the Young’s modulus as the only unknown. If the problem is discretized and formulated such that the Young’s moduli are the unknowns, with known displacements and strains, a forward solution technique can be applied to estimate the parameters [11–16]. Most of the methods in this category require the 2D displacement and strain fields as well as their first and second derivatives in the region of interest. However, due to the low lateral resolution of ultra- sound, the estimated lateral displacements have low signal-to-noise ratios (SNR). Once the low SNR of the lateral displacements has been overcome and the boundary forces have been measured, a reliable elasticity distribu- tion may be obtained [16]. In another approach, the inverse problem of elasticity is solved in the sense that a specific functional is minimized [5, 17–22]. Usually, this func- tional is considered to be the quadratic norm of the difference between the measured displacements and the displacements resulted from an assumed distribution of elasticity. Using the finite element method (FEM), itera- tive strategies based on Gauss-Newton or quasi-Newton methods have been proposed in the literature. Also, a method has been suggested to evaluate the gradient using adjoint equations in order to increase the computational efficiency [20]. The exact analytical solutions for the inverse problem of elasticity and other quasi-static parameters have been investigated in [23], where it was shown that the problem of determining the shear modulus be- comes unstable for nearly incompressible materials. It was shown that The uniqueness of the solution is shown to be highly dependent on the regulariza- 83 Chapter 4. Dynamic Finite Element Inverse Problem tion, available boundary conditions, prior knowledge of the elasticity on the boundaries and the incompressibility of the medium [24–26]. Recently, the authors have shown that applying wide-band excitation and calculating the asymptotic low-frequency magnitude of the transfer functions can increase the accuracy of the elasticity estimations [27]. Estimation of the viscosity, however, requires applying dynamic exci- tation to the material and calculating higher order derivatives (material velocities). Hence, compared with the elasticity reconstruction, the viscos- ity estimates will be more sensitive to noise. Viscoelastic modeling and parameter identification have also been explored extensively in the litera- ture [3, 6, 9, 10, 28–33]. Basically, through analysis of the attenuation and phase change of the propagating waves or the relaxation behavior of the structure, the dynamic properties can be characterized. The problem of re- constructing the damping matrix in a dynamic FEM framework has been mostly studied in the context of structural dynamics [34–36], where the goal is to identify an overall damping matrix for a given structure. In this work, the inverse problem of reconstructing the stiffness and viscosity using a dynamic finite element approach is tackled. A model is proposed that incorporates the Voigt model for the viscoelastic deformation of soft tissues in a dynamic finite element formulation. The inverse prob- lem of viscosity and elasticity is solved based on harmonic measurements of the axial displacements. The algorithm is devised such that by knowing the parameters on the boundaries, their distribution can be estimated in- side the medium, without needing to know the force on the boundaries or the lateral displacements in the region of interest. The method is similar to the elasticity reconstruction algorithm by Kallel and Bertrand in [17], in the sense that a displacement functional is minimized using a Gauss- Newton approach. However, the problem has been reformulated to account for the viscous damping effect, unknown lateral displacements and effects of the unknown forces on the Jacobian matrix. The accuracy of the pro- posed algorithm is evaluated using dynamic finite element simulations. Pre- liminary studies have been performed on tissue-mimicking phantoms with embedded elasticity and viscosity inclusions. Gelatin phantoms with con- trolled proportion of gelatin powder are widely used to model the mechanical and ultrasonic characteristics of human soft tissues. The estimation results are compared with other parameter identification techniques and rheometry data. 84 Chapter 4. Dynamic Finite Element Inverse Problem 4.2 Model 4.2.1 Linear Viscoelastic Model In this paper, normal fonts denote scalar parameters while vectors are shown in lowercase bold and matrices in uppercase bold. The overall stress tensor (σij) for a viscoelastic material can be divided into elastic and viscous terms as follows [37]: σij = σ elas ij + σ visc ij , (4.1) where, σelasij = λ ǫkkδij + 2µǫij , (4.2) σviscij = λ ′ǫ̇kkδij + 2µ ′ǫ̇ij , (4.3) Strain is denoted by ǫij and its time derivative by ǫ̇ij . Tensor notation is used, thus ǫkk = ǫ11+ǫ22+ǫ33 and δij = 1 if i = j and 0 if i 6= j. λ and µ are the Lamé constants and λ′ and µ′ are the viscosity characteristic parameters. In this paper, the viscosity constant of the medium is defined as η ∆ = 2µ′. The Lamé constants can be expressed in terms of the Young’s modulus (E) and the Poisson’s ratio (ν) of the medium. For the case of an isotropic material, all of the elastic parameters can be defined in terms of only two independent constants of the above. In most soft tissues, nearly static incompressibility is assumed where ν ≈ 0.5 and zero fluid compressibility results in λ′ = −η/3 [37]. With these assumptions, the matrix representation of (4.1) is as follows: σ = Cǫ+C′ǫ̇ , (4.4) where vectors σ and ǫ contain the six components of the isotropic stress and strain tensors and C and C′ can be derived from equations (4.1)-(4.3). C is the material elasticity characteristic matrix (see for example [37] or [38]) and C′ is the material viscosity characteristic matrix: C′ = η 3 2 −1 −1 0 0 0 −1 2 −1 0 0 0 −1 −1 2 0 0 0 0 0 0 3 0 0 0 0 0 0 3 0 0 0 0 0 0 3 . (4.5) 85 Chapter 4. Dynamic Finite Element Inverse Problem 4.2.2 Finite Element Model of the Deformation If the region of interest is discretized, by requiring the work of external forces to be equal to the work of inertia plus the work of the elastic and viscous forces, the dynamic finite element model of the deformation can be obtained. In this derivation, the internal deformation of the body, which is governed by equation (4.4), determines the work of the internal viscoelastic forces. The interested reader is referred to [39] for a detailed derivation. With the excitation f(t) and displacement u(t) for all the nodes as a function of time, the well-known transient FEM model can be obtained as follows: Ku(t) +Bu̇(t) +Mü(t) = f(t) , (4.6) where K is the stiffness matrix which depends on the elements elasticity vector e, B is the damping matrix and M is the mass matrix composed by globalization of the individual masses [39]. Note here that f is a vec- tor of the known boundary conditions, including forces and displacements, and therefore it depends on the elasticity distribution [40]. Usually M is diagonalized by assuming lumped masses at the nodes. If equation (4.5) is used in deriving (4.6), B will be a function of the viscosity vector η of the elements. In a general 3D problem with n nodes and m elements, K, B and M are 3n × 3n matrices, e and η are vectors of size m and u(t) and f(t) are vectors of size 3n. The Fourier transform of equation (4.6) yields:( K + jωB − ω2M) û = f̂ . (4.7) where hatted variables denote Fourier transforms and j is the imaginary unit. Decomposing û and f̂ into their real and imaginary parts, equation (4.7) can be written as: A [ ûr ûi ] = [ f̂ r f̂ i ] , (4.8) with A ∆ = [ K − ω2M −ωB ωB K − ω2M ] , (4.9) where the superscripts r and i denote the real and imaginary components of a vector, respectively. 86 Chapter 4. Dynamic Finite Element Inverse Problem 4.3 Inverse Problem Upon applying an external excitation to a soft material, a deformation will be produced within the material, which depends on the mechanical proper- ties, boundary conditions and the geometry. A method was proposed in [17] to optimize for the tissue Young’s modulus in a finite element model, such that the static least squares error between measured and model-predicted displacements is minimized. Similar optimization procedures have been studied in a number of other works with static or quasi-static displace- ment data [5, 18–22]. The idea is to minimize a functional that is based on the difference between the observed real-valued displacement vector and the predicted one. If, however, a transient or harmonic excitation is applied, the motion will also be affected by dynamic properties according to (4.6) or (4.7). Therefore, a more general formulation needs to be derived that takes into account the real and imaginary displacement components and optimizes for the best viscoelastic parameters that match the model. 4.3.1 Reconstruction Method An inverse algorithm is devised that searches for the optimal values of the Young’s modulus and viscosity that yield the least deviation from the ob- served displacements inside the medium. With measured displacement vec- tor û0, an error functional can be defined as follows: φ(p) = 1 2 ‖ûr(p)− ûr0‖2 + 1 2 ∥∥ûi(p)− ûi0∥∥2 , (4.10) where p = [eT ηT ]T is a 2m parameter vector consisting of the Young’s moduli and viscosities of the m elements in the model. The parameter p can be iteratively updated using the Gauss-Newton algorithm or any descent method. In this paper, the Jacobian has been calculated by differentiating equation (4.8) and the step direction has been computed with the Levenberg-Marqardt algorithm. The parameter vector is updated at each iteration as pk+1 = pk +∆pk, k ∈ N, where k is the iteration index. The update vector, ∆pk, can be calculated from: Hk∆pk = −JkT∆ûk , (4.11) where Jk is the Jacobian, Hk is the Hessian and ∆ûk is defined as follows: ∆ûk = [ ∆ûrk ∆ûik ] = [ ûrk − ûr0 ûik − ûi0 ] , (4.12) 87 Chapter 4. Dynamic Finite Element Inverse Problem where, for simplicity, ûk denotes û(pk). Note that for a 3D FEM mesh with n nodes, ∆ûk is a 6n parameter vector. 4.3.2 Calculation of the Jacobian and the Hessian The Jacobian matrix describes the sensitivity of the displacements with respect to the parameters, and is defined as follows: J ∆ = ∂û r ∂p ∂û i ∂p = ∂û r ∂e ∂û r ∂η ∂û i ∂e ∂û i ∂η . (4.13) Differentiating (4.8) with respect to the jth element of p, i.e. pj , one obtains: ∂A ∂pj [ ûr ûi ] +A{J j} = ∂f̂ r ∂pj ∂f̂ i ∂pj , (4.14) where {J j} denotes the jth column of J defined in (4.13). Note that pj can be either the elasticity or the viscosity of an element, therefore ∂A/∂pj can be computed from (4.9) as follows: ∂A ∂pj = [ ∂K ∂ej 0 0 ∂K ∂ej ] , pj = ej ω [ 0 −∂B ∂ηj ∂B ∂ηj 0 ] , pj = ηj (4.15) The right-hand side of equation (4.14) can be calculated by considering the effect of the displacement boundary conditions on the vector f̂ . As a result, with knownA, equation (4.14) can be solved for {J j} and thus the Jacobian matrix can be constructed. Once the Jacobian has been evaluated at the kth iteration, first-order approximation of the Hessian yields: Hk = Jk TJk . (4.16) Note that this modified Hessian, ignores the effect of the second order resid- uals in the least squares problem of (4.10). However, in the proximity of the 88 Chapter 4. Dynamic Finite Element Inverse Problem solution, this error is negligible [41]. In general, with n nodes and m elements with unknown parameters in the mesh, J is a 6n× 2m matrix and H is a 2m× 2m matrix. 4.3.3 Practical Considerations Equation (4.11) requires the Hessian to be invertible. However, when the Jacobian is rank deficient, the Hk in (4.16) will be badly conditioned and may not be easy to invert. The Levenberg-Marquardt algorithm modifies the Hessian in such a way that the resulting step will be a combination of the directions predicted by the steepest descent and the Gauss-Newton methods, depending on how far the current point is from a local minimum. If far from a local minimum, the algorithm behaves like the steepest descent method. However, in the vicinity of a minimum, the step direction will be close to that of the Gauss-Newton method and the positive definiteness of the Hessian can be gauranteed [41, 42]. The Hessian may thus be modified as follows: Hk = J T k Jk + λkI2m×2m , (4.17) where I2m×2m is an identity matrix with m being the number of elements with unknown parameters. λk is a small regularization factor that makesHk positive definite. This method has been reportedly used with success in the context of the inverse problem of elasticity (e.g. [17]). The required update vector can be obtained by utilizing a trust-region or line-search procedure [41]. The reconstruction result is highly dependent on the displacement noise. As suggested by Doyley et. al a spatial filter can be applied on the modulus distribution at every iteration to make the solution of the problem smooth [43]. As a result, the modulus of each element will be a weighted sum of its own value and the moduli of its adjacent elements. A linear filter can thus be constructed in the form of a sparse matrix that contains the required weights for each element. The filter can be convolved with the updated modulus distribution at each iteration to conduct the optimization toward a smooth solution. To constrain the problem and avoid physically infeasible or implausible results, assumptions can be made on the parameters by adding constraints to (4.10). In particular, negative elasticity and viscosity values should be avoided. To bound a parameter p to values greater than or equal to pL, a technique based on the change of variables can be implemented. Therefore, an auxiliary variable s can be defined such that p = s2+pL. Using the chain 89 Chapter 4. Dynamic Finite Element Inverse Problem rule, the problem can be formulated and solved for s, and finally p can be calculated from the optimized value of s. To ensure a sufficient decrease of the functional at every iteration and guarantee convergence, a line search procedure has been implemented to determine the step-size based on the Armijo condition which uses a first order approximation of the functional [41]. 4.4 Tests 4.4.1 Two Dimensional Modeling Plane-Stress vs. Plane-Strain While a 3D model accounts for the deformations in all directions, certain assumptions can be made to simplify the problem to only two dimensions. When the medium is confined such that the out-of-plane strain and displace- ments are minimized, a plane-strain model is used to analyze the deforma- tions. This happens when some walls parallel to the imaging plane confine the medium or when the thickness of the medium is significantly larger in the elevational direction. Also if the axial compression is applied uniformly on the top surface of the medium, a plane-strain assumption may be valid on the middle plane [7]. On the other hand, a plane-stress assumption can be legitimately made when the elevational extent of the medium is much smaller than the other two dimensions. In this work, for error quantification and performance analysis, a 2D plane-stress case is assumed. Under this assumption, the material charac- teristic matrices for the elasticity and viscosity can be computed by setting the out-of-plane components of the elastic and viscous stress tensors equal to zero. This results in the following 2D characteristic matrices: C = E (1− ν2) 1 ν 0ν 1 0 0 0 (1−ν)2 , (4.18) C′ = η 2 1 −1 0−1 1 0 0 0 2 , (4.19) where E, η and ν are the Young’s modulus, viscosity and Poisson’s ratio, re- spectively. Alternatively, under plane-strain assumption, the characteristic 90 Chapter 4. Dynamic Finite Element Inverse Problem Figure 4.1: 2D finite elements model, using rectangular elements. The region is bounded at the top and excited from below by a force or displacement distribution of a known frequency. Note that the arrows do not specify the nodes that are excited, but the presence of a distributed or localized excitation. The axial and lateral directions of the ultrasound are shown. 91 Chapter 4. Dynamic Finite Element Inverse Problem Assumean initial distribution for p Run forward problem at frequency ù Calculate axial displacements û (p) and the Jacobian a For simulated model or experiment, measure axial displacements û at0a ù Using the Levenberg-Marqardt method, update the parameters p Calculate the functional and its gradient using only axial displacements ö(p) puJ ¶¶= /ˆ a Convergence? Finish Yes No Figure 4.2: Block diagram of the optimization procedure. 92 Chapter 4. Dynamic Finite Element Inverse Problem Width (mm) D ep th (m m) Locations of the Inclusions 0 10 20 30 40 0 10 20 30 40 Elasticity Inclusion Viscosity Inclusion (a) Width (mm) D ep th (m m) Real Part of Strain 0 10 20 30 40 0 10 20 30 40 (b) Width (mm) D ep th (m m) Imaginary Part of Strain 0 10 20 30 40 0 10 20 30 40 (c) Width (mm) D ep th (m m) Real Part of Strain 0 10 20 30 40 0 10 20 30 40 (d) Width (mm) D ep th (m m) Imaginary Part of Strain 0 10 20 30 40 0 10 20 30 40 (e) Width (mm) D ep th (m m) Estimated Elasticity (kPa) 0 20 40 0 10 20 30 40 15 20 25 (f) Width (mm) D ep th (m m) Estimated Viscosity (Pas) 0 20 40 0 10 20 30 40 15 20 25 (g) Figure 4.3: The axial strain around an elasticity and a viscosity inclusion. The locations of the inclusions are shown in (a). The real and imaginary parts of axial strain due to a 5Hz harmonic displacement excitation of the lower surface are in (b) and (c). The real and imaginary axial strain profiles due to applying the same excitation to only a few nodes at the middle bottom of the phantom are in (d) and (e). A 2D dynamic FEM model has been used to generate the strain fields. The artifacts in the strain images are caused by several factors such as the loading profile, boundary conditions and the geometry. By solving the inverse problem, the elasticity and viscosity inclusions can be accurately delineated as illustrated in (f) and (g). The inverse problem was solved for the case in (b) and (c). Similar results are obtained with the data from the other case. 93 Chapter 4. Dynamic Finite Element Inverse Problem matrices change to the following ones: C = E (1 + ν)(1− 2ν) (1− ν) ν 0ν (1− ν) 0 0 0 (1− 2ν) , (4.20) C′ = η 3 2 −1 0−1 2 0 0 0 3 . (4.21) Using Only Axial Displacements If the lateral displacements are not available, the functional in equation (4.10) can be evaluated in terms of the axial data only. Therefore, once the Jacobian has been calculated, the rows that correspond to the lateral data can be removed and then the step can be evaluated from equation (4.11). In the following simulations and experimental analysis, only axial displace- ments were used in the inverse problem and the elasticity and viscosity at the top of the phantom were assumed constant, thus constraining the op- timization. The lateral displacement and force distribution were not used in solving the inverse problem. The block diagram for the procedure to solve the inverse problem in the simulation and experiments is depicted in Fig. 4.2. A stopping criterion monitoring the gradient size and the number of iterations has been used. Three criteria were checked in order to stop the iterations, including the value of the functional, the functional gradient and maximum number of iterations. The thresholds for these conditions were chosen by trial and error. 4.4.2 Numerical Simulations A 2D 4cm × 4cm region, bounded from above with laterally nonslip con- ditions, has been meshed with four-node rectangular elements. Models of the elasticity and viscosity as defined in equations (4.18) and (4.19) were adopted. Figure 4.1 shows a typical finite element grid as used in the simu- lations. As can be seen, a rectangular region of interest is bounded at the top and excited from below by applying force or displacements to one or several nodes. In the plane-stress formulation, the thickness of the region has been assumed to be 1cm. Subject to a force or displacement excitation, the nodal deformation could be calculated as explained in Section 4.2.1. With a static displacement excitation, it has been shown that the solution to the inverse problem of elasticity in the plane-strain state is unique, given the value of 94 Chapter 4. Dynamic Finite Element Inverse Problem the Young’s moduli on the lower or upper boundaries of the medium [24]. The uniqueness of the inverse problem of viscoelasticity under plane-stress dynamic deformation has not been studied explicitly in the literature. It is possible that the same plane-strain analysis of [24] can be performed for the plane-stress case to obtain similar uniqueness criteria with complex-valued parameters as in our case. We proceed with the assumption that a unique solution can be obtained by knowing the viscosity and elasticity on the top row of the finite element grid, given only the complex axial displacements. A constant density of 1000kg/m3 is assumed for the medium which is typical of soft tissues. Within constant regions, the Young’s modulus and viscosity are assumed to be 10kPa and 10Pas, respectively, which are chosen arbitrarily in the range of the viscoelastic properties of human soft tissues. The Poisson’s ratio has been assumed to be 0.495 to mimic the near incom- pressibility of soft tissues. To elucidate the significance of modulus imaging versus interpretation of the strain images, a medium with an elasticity and a viscosity inclusion has been simulated. Figure 4.3(a) shows the locations of the two inclusions. The background has a Young’s modulus of 10kPa and a viscosity of 10Pas, while in the elasticity inclusion the Young’s modulus is 30kPa and the vis- cosity of the viscous inclusion was 30Pas. A symmetric mesh of 21 × 21 nodes has been applied to the region and a 5Hz steady-state compressional displacement excitation has been applied to the phantom in the axial direc- tion. In a first trial, all of the bottom nodes of the phantom were excited and the resulting axial strains were computed. Figs. 4.3(b) and 4.3(c) show the real and imaginary parts of the axial strain, respectively. In a second test, only five middle nodes at the bottom of the grid were excited, which resulted in the strain images in Fig. 4.3(d) and 4.3(e). One can notice a high dependency of the strain images on the loading and boundary condi- tions. The presence of artifacts in the images makes it hard to delineate the inclusions by direct interpretation of the strain images, however, solving an inverse problem with the knowledge of the boundary conditions and the loading profile enables one to reconstruct the modulus distributions accu- rately. Using the displacements under the first described loading conditions, the reconstructed images are shown in Fig. 4.3(f) and 4.3(g) for the Young’s modulus and viscosity. Due to zero noise in this simulation, the parameters did not need smoothing, therefore perfect reconstruction can be observed. Similar images could be obtained using the data from the second loading condition. The sensitivity of the estimations to the displacement noise is illustrated in Fig. 4.4. The same material properties and inclusion locations as in 95 Chapter 4. Dynamic Finite Element Inverse Problem 20 40 60 80 100 10−4 10−3 10−2 10−1 Displacement SNR (dB) Es tim at io n Er ro r Elasticity RMS Error Not Filtered Filtered (a) 20 40 60 80 100 10−4 10−3 10−2 10−1 Displacement SNR (dB) Es tim at io n Er ro r Viscosity RMS Error Not Filtered Filtered (b) 20 40 60 80 100 10−1 100 101 102 103 Displacement SNR (dB) Es tim at io n CN R Elasticity CNR Not Filtered Filtered (c) 20 40 60 80 100 10−1 100 101 102 103 Displacement SNR (dB) Es tim at io n CN R Viscosity CNR Not Filtered Filtered (d) Figure 4.4: The displacement signal-to-noise ratio (SNR) changes the accu- racy of the estimation. The computed displacements have been corrupted by a known power of noise and the parameters have been estimated using the raw data (the results are in black solid lines) or spatially filtered dis- placements (the results are in gray dashed lines). The RMS error for the elasticity and viscosity estimates are shown in (a) and (b) respectively, while the contrast-to-noise ratio (CNR) for those estimates are depicted in (c) and (d). 96 Chapter 4. Dynamic Finite Element Inverse Problem Fig. 4.3(a) have been used and the same mesh with 21× 21 nodes has been applied to the region. The entire lower side of the region was subjected to a 5Hz steady-state excitation. White noise was then added to the resulting displacements and the inverse problem was solved for the elasticity and viscosity values. Two criteria were utilized to assess the performance of the estimation. First, the RMS error was measured based on the difference between estimated moduli and the actual ones. The RMS error for the elasticity can be defined as follows: ǫ(e) = √√√√ 1 m m∑ i=1 ( eei − e0i )2 (e0i ) 2 , (4.22) where m is the number of the elements and e0i and e e i are the actual and estimated elasticities for the ith element. The RMS error for the viscosity is defined in the same way. For each level of SNR, the model was simulated 50 times and the means and standard deviations of the RMS error for elasticity and viscosity are shown in Fig. 4.4(a) and 4.4(b). If instead of using the raw noisy displacements, a small 3× 3 averaging filter is applied to the displace- ments prior to running the inverse algorithm, a significant improvement can be achieved in lower displacement SNR values (around 40dB), while at high SNRs, the smoothing results in poor estimation at the boundaries of the inclusion, and thus a higher RMS error. As a second criterion to analyze the exactness of the solution, the contrast- to-noise ratio (CNR) of the estimated elasticity and viscosity images with inclusions can be used [45]. If the area of interest is composed of two homo- geneous regions with different Young’s moduli, the CNR(e) can be defined as follows: CNR(e) = 2(s1 − s2)2 σ21 + σ 2 2 , (4.23) where s1 and s2 are the average elasticities estimated at the inclusion and at the background and σ1 and σ2 are the standard deviations of the elasticity estimates in those regions. CNR(η) can also be defined for the viscosity estimates in a similar way. A high CNR value will be obtained if the esti- mates are smooth in each region but quite discriminated between the two. The means and standard deviations of the CNR values versus displacement SNR for the aforementioned phantom are plotted in Fig. 4.4(c) and 4.4(d) for the elasticity and viscosity. It can be seen that the contrast is enhanced when the displacements are filtered prior to solving the inverse problem at higher noise levels, while the contrast is compromised at higher SNR values. 97 Chapter 4. Dynamic Finite Element Inverse Problem Width (mm) D ep th (m m) Problem Geometry 0 20 40 60 0 10 20 30 40 50 60 Elasticity Inclusion Viscosity Inclusion Region of Interest Figure 4.5: Solving the inverse problem for a region of interest (ROI) with unknown lateral boundary conditions. A 6cm× 6cm region has been simu- lated while only the axial displacement at a 4cm×4cm region in the middle has been used for parameter identification. 98 Chapter 4. Dynamic Finite Element Inverse Problem Width (mm) D ep th (m m) Elasticity (kPa) 0 10 20 30 0 10 20 30 10 15 20 25 30 (a) Width (mm) D ep th (m m) Elasticity (kPa) 0 10 20 30 0 10 20 30 10 15 20 25 30 35 40 (b) Width (mm) D ep th (m m) Elasticity (kPa) 0 10 20 30 0 10 20 30 10 15 20 25 30 (c) Width (mm) D ep th (m m) Viscosity (Pas) 0 10 20 30 0 10 20 30 10 15 20 25 30 (d) Width (mm) D ep th (m m) Viscosity (Pas) 0 10 20 30 0 10 20 30 0 10 20 30 (e) Width (mm) D ep th (m m) Viscosity (Pas) 0 10 20 30 0 10 20 30 0 10 20 30 40 (f) Figure 4.6: The effect of unknown boundary conditions in the reconstruction of elasticity and viscosity. The medium in Fig. 4.5 has been simulated. The axial displacements were used to solve the inverse problem in the illustrated ROI. Different lateral boundary conditions were assumed and the inverse problem was solved: no lateral boundary forces were assumed in (a) and (d); the top and the bottom of the ROI were assumed not to move laterally in (b) and (e); and all four sides of the ROI was assumed to have zero lateral motion in (c) and (f). Note that the simulations were performed for an FEM model with 21× 21 elements, while the displayed images are upsampled by a factor of 5 for better illustration. 99 Chapter 4. Dynamic Finite Element Inverse Problem If the lateral displacements and forces at some or all parts of the bound- aries are not available, some artifacts may be present in the reconstructed parameters. Figure 4.6 illustrates the effect of unknown boundary conditions on the reconstruction results. The same medium as depicted in Fig. 4.3(a) has been embedded in the middle of a larger region with the same back- ground properties. The forward problem has been solved for a 6cm × 6cm region. As in the previous case, the lower end of the phantom was excited by a 5Hz axial harmonic motion, while the top was held stationary. The region of interest (ROI) was the middle 4cm × 4cm of the phantom that enclosed the two inclusions as shown in Fig. 4.5. The parameters have been recon- structed using only the axial displacements in the ROI and with different assumptions on the lateral motion at the boundaries. Spatial filtering has not been used on the estimated parameters in order to obtain a clear com- parison. First, zero lateral force has been assumed on the boundaries. The reconstructed images of the elasticity and viscosity are shown in Fig. 4.6(a) and 4.6(d). Next, the top and the bottom of the region have been assumed to have zero lateral motion, while the right and the left boundaries were free to move laterally. The results are shown in Fig. 4.6(b) and 4.6(e). Finally, all four sides of the region were assumed to have zero lateral motion and the results are depicted in Fig. 4.6(c) and 4.6(f). Another source of error in solving the inverse problem of viscoelasticity arises from the fact that the deformation of the body is three dimensional. Plane-stress or plane-strain assumptions are intended to simplify the model to 2D, however such simplifications may result in inaccurate estimations. To demonstrate this, a 3D region (4cm× 4cm× 3cm) has been simulated with the same properties and inclusions as denoted before. Having the top of the phantom axially fixed and free to move in the lateral direction, a 10Hz displacement excitation has been applied to the bottom surface in the axial direction. The forward problem was solved using equation (4.7). As illus- trated in Fig. 4.7(a), the axial displacements in the middle plane were used to reconstruct the elasticity and viscosity using the proposed method. The re- construction was performed once with plane-stress assumption (Figs. 4.7(b) and 4.7(d)) and another time using plane-strain matrices (Figs. 4.7(c) and 4.7(e)). The results show that a 2D plane-stress approximation can be used to model such a geometry, however the artifacts can be significantly reduced if a plane-strain assumption is applied. Besides the 2D simplification of 3D deformations, a significant error may arise from using an inaccurate model in the procedure. To show how the reconstruction result may be affected if a different structure of the damping matrix is assumed, the test with a 5Hz distributed displacement excitation 100 Chapter 4. Dynamic Finite Element Inverse Problem y x z Regionof Interest (Analysis Plane) Simulated 3D Phantom Elasticity Inclusion Viscosity Inclusion (a) Width (mm) D ep th (m m) Estimated Elasticity (kPa) 10 20 30 40 10 20 30 40 10 15 20 (b) Width (mm) D ep th (m m) Estimated Elasticity (kPa) 10 20 30 40 10 20 30 40 10 15 20 25 30 (c) Width (mm) D ep th (m m) Estimated Viscosity (Pas) 10 20 30 40 10 20 30 40 10 20 30 40 50 60 (d) plane-stress Width (mm) D ep th (m m) Estimated Viscosity (Pas) 10 20 30 40 10 20 30 40 10 20 30 40 50 (e) plane-strain Figure 4.7: The effect of the 3D deformations on the estimation when a 2D model is used to reconstruct the parameters. (a) shows the 3D phantom that has been simulated with a high elasticity and a high viscosity inclusion. The axial displacements in the middle plane were used to solve the inverse problem in 2D. The reconstructed elasticity and viscosity with plane-stress assumption are shown in (b) and (d), while having a plane-strain assumption produces elasticity and viscosity estimates in (c) and (e), respectively. 101 Chapter 4. Dynamic Finite Element Inverse Problem on the bottom of the phantom has been repeated using the same geom- etry and inclusions as before. To simulate the displacements at the first stage, the damping matrix has been constructed in a similar way to the stiffness matrix, i.e. equation (4.18) was used instead of (4.19), with a mul- tiple of the viscosity substituting E0. This is in accordance with Rayleigh damping where the damping matrix of every element is proportional to its stiffness matrix [46]. This model of damping, however, has not been taken into account in the reconstruction, thus the original damping structure as it results from equation (4.19) has been used in the inverse procedure. The estimates as shown in Fig. 4.8 indicate that the inclusions can still be iden- tified; however, the presence of artifacts in the viscosity image are noted as a consequence of improper model assumptions. Some correlation can also be noticed between the artifacts in the viscosity image and the elasticity distribution. 4.4.3 Experiments A tissue mimicking phantom has been constructed with 12% (by weight) bovine skin gelatin in water. 2% cellulose (by weight) has been uniformly added to the phantom prior to solidification to act as ultrasound scatterers and then the mixture was poured in a 38(x)×40(y)×25(z) mm3 mold. A small piece of a polyvinyl alcohol (PVA) sponge (Ceiba Technologies, Chan- dler, AZ, USA) was embedded in the phantom shortly before solidification. The sponge was soaked in water and degassed. This type of PVA sponge has a relatively high Young’s modulus compared to its background gelatin and a significantly higher relaxation-time and viscosity [8]. The phantom was placed on a specially designed shaker [8] and ultrasound RF data were captured. The axial and lateral motions of the phantom at the transducer side were nearly zero while the axial displacement of the other side was estimated. Since the phantom was ultrasonically imaged up to a depth of 37mm (before the bottom of the phantom), the lateral motion at the bot- tom of the phantom is not known. It was assumed that lateral force on the phantom was insignificant. A wide-band displacement excitation (1–30Hz) with a Gaussian amplitude distribution with standard deviation of 44µm was applied to the bottom of the phantom. The internal motion of the phantom was estimated using a cross-correlation based algorithm with prior estimates [47]. Frequency analysis was performed on the displacements and the real and imaginary parts at 10.7Hz were selected to solve the inverse problem. Since power spectral analysis yields spectral data at discrete fre- quencies depending on the window length that is applied on the time-domain 102 Chapter 4. Dynamic Finite Element Inverse Problem data [48], only data at certain frequencies were available to solve the inverse problem. Therefore, displacements at 10.7Hz were chosen arbitrarily in the frequency range of excitation, such that considerable viscoelastic behavior could be observed in the phantom. In conventional ultrasound, an image is acquired by collecting data from every line, progressively. Hence, there is a systematic delay between the data recorded from different lines of the image, which causes a linear phase difference between the displacements at different lines. As a pre-processing stage on the experimental data, this linear lateral phase gradient due to the limitation of the acquisition time between different lines were calculated and its effect was removed. For this reason, the phase of the displacements were measured and a linearly increasing trend was optimally identified in their lateral profiles. By removing this trend, the acquisition delays were com- pensated. The displacements were originally estimated in a grid of 49 lines by 67 blocks, filtered and then down-sampled by a factor of 2. As shown in Fig. 4.4, filtering the displacements improves the accuracy of the estima- tion at the range of SNR that is typical for elastographic motion estimation. The down-sampling was performed to reduce the number of variables and the size of the problem. This has only achieved through compromising spa- tial resolution with processing speed and system memory. Based on independent rheometric measurements, the background gelatin had a Young’s modulus and a relaxation-time of approximately 15kPa and 1ms respectively. The relaxation-time is the ratio of the viscosity to the Young’s modulus of the material. The relaxation-time of the PVA sponge was also measured to be approximately 4.5ms at 10.7Hz [8]. Using initial conditions of 15kPa and 15Pas for the elasticity and viscosity for the entire medium, the inverse problem was solved and the parameters were recon- structed. Since the elevational size of the phantom was nearly one half of the other two dimensions, a plane-stress model was presumably more accu- rate than a plane-strain one in describing the internal deformations. The iterative scheme was set to stop after the gradient reached a sufficiently small value [41]. The B-mode image and approximate location of the inclusion are depicted in Fig. 4.9(a). The reconstructed images of elasticity and viscosity are shown in Fig. 4.9(b) and 4.9(c). The estimated relaxation-time is also shown in Fig. 4.9(d). It can be seen from the reconstruction results that the relaxation-time differ- ence between the medium and the inclusion is approximately 2ms. This parameter was measured to be 3.5ms in rheometry tests. Although some artifacts are present in the estimated images, the inclusion can be clearly distinguished from the background material. 103 Chapter 4. Dynamic Finite Element Inverse Problem 4.5 Discussion The feasibility of reconstructing the elasticity and viscosity in a non-homogeneous medium based on a finite element analysis has been studied in this work. A model has been devised to account for the viscoelastic changes in soft tissues. The simulations and experiments were performed in 2D while they can be generalized to a 3D problem if other components of the displacements can also be measured experimentally. The assumption of a nearly plane-stress deformation has been made in this work and the simulations and inversion techniques have been utilized accordingly. If, however, a plane-strain state is known to explain the observed deformations more accurately, the same inverse algorithm can be used with a proper plane-strain formulation of the forward problem. While 2D FEM is a simplification of the more general 3D problem, even the proposed 3D formulation or other 3D models in the literature may not accurately model the deformations. As explained in Section 4.2.2, the mass matrix in dynamic finite element analysis can be diagonalized to improve computational efficiency. A diago- nal representation of the damping matrix is also permissible in some situa- tions. A more general scheme to characterize the damping matrix, coined as Rayleigh or proportional damping, is to form it as a linear combination of the stiffness and mass matrices [46]. Depending on the structure and frequency, the effect of either stiffness or mass matrices may be dominant, which re- spectively make the damping matrix approximate a consistent or lumped representation. The damping matrix using this approach is frequency de- pendent, where the part attributable to the stiffness matrix increases with increasing frequency and the part attributable to the mass matrix increases with decreasing frequency [39]. A more realistic characterization of the damping matrix for soft viscoelastic materials has been proposed in this paper. The equations of dynamic equilibrium have been discretized in a finite element mesh and further simplified to obtain the desired relationship between the force and displacement vectors. With this representation, one can be sure that within the discretization error, the acquired finite element model agrees with other viscoelastic formulations based on the equilibrium equations. To reduce the condition number of the Hessian matrix (often larger than 1012), two methods have been tested. In the Levenberg-Marquardt method, the step is close to that of the Gauss-Newton when the Hessian is invertible and similar to that of the steepest descent when Hessian is non-invertible. Another technique is to modify the small eigenvalues of the problem to limit the condition number of the Hessian with an upper bound. In the 104 Chapter 4. Dynamic Finite Element Inverse Problem simulations and experiments in this paper, the difference between the results from the two methods was seen to be insignificant, thus the Levenberg- Marquardt technique has been selected. The poor condition numbers of the inverse problems in elasticity and viscosity make the final solution dependent on the initial conditions and the optimization technique that is used. One approach to obtain smooth profile of the parameters is to add the Laplacian of the parameters to the functional and search for the solution that minimizes the displacement difference and the Laplacian term at the same time [20]. In this paper, as suggested in [43] the parameters were spatially filtered by a 3 × 3 averaging kernel at each iteration to render a smooth solution when dealing with noisy displacements. The inverse problem of elasticity and viscosity is solved using only axial displacements. Figure 4.3 demonstrates the artifacts in the strain images and how they can be removed by solving the inverse problem. While knowl- edge of the lateral boundary conditions will result in accurate estimation and removal of the artifacts, from Fig. 4.6 it can be seen that axial data are not sufficient to characterize the medium. With lateral motion estimation in elastography, axial and lateral displacements can both be used to remove the artifacts corresponding to unknown boundary conditions. The simulations were performed with a relatively coarse mesh of 21× 21 within a 40×40mm2 region. The spatial resolution can be enhanced by using a finer grid size, at the expense of computational resources. Alternatively, spatial resolution can be enhanced recursively through mesh refining. For this purpose, a coarse mesh can be first applied to the original problem. Once the parameters have been reconstructed, the resolution can be enhanced in different regions of interest by applying a finer localized mesh, so that the inverse problem can be solved with higher spatial resolution only in that region. Based on Figs. 4.5 and 4.6, such a refining window can be applied to any part of the image, given that the lateral boundary conditions are partially known for that window. The error in estimation has been measured by defining two criteria, RMS error and contrast-to-noise ratio. Figure 4.4 shows that below a certain level of displacement SNR, parameter reconstructions results in noisy and unre- liable estimates. This sensitivity to noise may be lowered by appropriate filtering of the noisy displacement data. The RMS error seems to reach a plateau rather than converging to zero at high SNR values. This is due to filtering of the parameters at every iteration as explained in Section 4.3.3. Based on equation (4.22), the RMS error may originate from either estima- tion bias or variance [49,50]. The errors quantified in Fig. 4.4 for simulated data are mainly due to the latter, which indicates the lack of precision in the 105 Chapter 4. Dynamic Finite Element Inverse Problem estimates. Compared to the estimation variance, the bias can be made small by choosing suitable initial conditions for the recursions and assuming cor- rect parameters as the constrained boundary values for the problem. Both of these errors are important factors in producing reliable images; however, the presence of a small amount of bias in the estimated parameters may still yield decipherable images while estimation variance may drastically drown out the modulus distribution. In in vivo situations, it may not be feasible to know the exact boundary conditions of the tissue. In that case, lateral displacement measurements would help reduce the artifacts in the reconstructed parameters. However, due to the coarse lateral resolution of ultrasound imaging, if accurate lat- eral motion estimation is not obtainable, some basic assumptions should be made on the lateral force or displacement on the boundaries of the imaging window. Figure 4.6 shows how such assumptions influence the reconstruc- tions. In this example, uncertainty of the lateral boundary conditions had negligible effect on the estimated Young’s modulus, while estimated viscos- ity is seen to be dependent on the boundary conditions assumed. In this particular case, one can see that, in comparison to the other assumptions, a zero force assumption produced a more acceptable viscosity estimate in Fig. 4.6(d) while a faded shadow of the elasticity inclusion is visible. Al- though accurate knowledge of the lateral boundary conditions is necessary, it is not sufficient to produce estimations that are free of artifacts. As seen in Fig. 4.7 simplifying a 3D problem to a 2D one can also produce errors. In the case of Fig. 4.7(a), the elevational dimension of the phantom was comparable to the other dimensions, hence a plane-strain model was rela- tively more accurate in describing the deformations. Here, a plane-strain model would be best suited if the thickness of the phantom were increased, while a smaller thickness would make the result of a plane-stress model more reliable. The systematic estimation error caused by choosing inaccurate models are depicted in Figs. 4.7 and 4.8. Ideally, the reconstruction results in both of these cases should be the same as those in Figs. 4.3(f) and 4.3(g). Us- ing a 2D model to describe the 3D deformations will be acceptable only if the dimensionality reduction to obtain a plane-stress or plane-strain state is reasonable. The blurs and artifacts in Figs. 4.7(b) and 4.7(d) are due to the marginal accuracy of the plane-stress approximation. In this case, a plane-strain model produced less artifact and more clear images as seen in Figs. 4.7(b) and 4.7(d). Also based on the findings in Fig. 4.8, one should expect to see some artifacts from misrepresentation of the damping phenom- enon, since the model proposed by incorporating equation (4.5) or (4.19) into 106 Chapter 4. Dynamic Finite Element Inverse Problem the FEM equations may not accurately describe the damping mechanisms in soft tissues. The algorithms have been tested with experimental data. A gelatin- based phantom was constructed with a hard and more viscous sponge-like PVA inclusion. The estimation results are in agreement with the rheometry tests reported in [8]. However, rheological models do not account for tissue compressibility and compressibility was not considered in the FEM model either. During phantom construction, the small piece of sponge was degassed and saturated with water, so that the volumetric changes are likely to be dominated by water. Since the volume within the inclusion was preserved with the water trapped inside, a nearly incompressible condition has been assumed. To what extent this was achieved is an acknowledged challenge. The porosity of the sponge material may make it compressible and this would introduce errors in the parameter identification. However, in such an FEM model, it is only the viscoelasticity that can account for the observed phase changes, while porosity and compressibility result in a change in the static deformation profile. Lateral forces were not measured at the boundaries and maintaining non-slip lateral conditions at the boundaries is not trivial. Therefore, un- known lateral boundary forces and displacements in the problem produced artifacts in the reconstructed images as expected from the simulations. The gray shadows below the inclusion in both images of Fig. 4.9 are due to the lack of information about the boundary conditions as well as inaccuracy in estimating the displacement phasors. Due to lower sonographic SNR in the region below the sponge, relatively poorer displacement estimates were ob- tained, which also lowered the estimation accuracy in that region compared to the rest of the phantom. If lateral displacements are measured or lateral forces are known at the boundaries, the accurate information can be used to enhance the reconstruction results. In addition, a plane-stress assumption has been made in this experiment due to the relatively smaller elevational dimension of the phantom compared to the other dimensions. While plane- stress should be a more suitable model than plane-strain, it would be more accurate if the thickness were smaller. Therefore, some systematic error should be expected due to utilizing a simplified model. In in vivo experiments where boundary conditions are not known, high quality estimates of the axial and lateral displacements are necessary to solve the inverse problem in any region of interest. Given the much more complex environment when dealing with human subjects, knowledge of the lateral displacements may be crucial to decouple the effect of the boundary conditions. However, for simple situations where legitimate assumptions can 107 Chapter 4. Dynamic Finite Element Inverse Problem be made for the lateral displacements, the axial component of the motion may be sufficient to reconstruct the viscoelastic properties. 4.6 Conclusion The feasibility of measuring the viscoelastic parameters of soft tissue using dynamic finite element models has been explored in this paper. A finite el- ement model has been derived to account for the viscoelastic changes in the soft tissue, consistent with the known rheological models. Using this defor- mation model and the observed dynamic data, the viscosity and elasticity have been reconstructed in the medium by solving an inverse problem. The proposed algorithm was derived for the general 3D finite element model. The 2D algorithms were also explored to comply with conventional ultrasound elastography. It is known that the solution to the elasticity inverse problem using quasi-static FEM is subject to high noise sensitivity. The accuracy of estimating the viscoelastic parameters using the proposed algorithm was found to be dependent on the boundary conditions, displacement noise level and appropriate temporal and spatial filtering of the displacements. An it- erative solution to the inverse problem of elasticity and viscosity has been formulated. A well-known Gauss-Newton based approach to solve the in- verse problem of elasticity has been modified in such a way that steady-state harmonic deformations are used to estimate the viscoelastic properties and the need for lateral information is minimized. The solution can be computed through incorporating a priori knowledge of the moduli, such as the values of the parameters on the boundaries, the smoothness of the solution, etc. The methods were studied in FEM simulations and the effect of the displace- ment noise and boundary conditions were explored. The approach has been tested in an experiment and successfully reconstructed the viscosity and elasticity in a gelatin based phantom with an inclusion of known properties. It was found that in some simple situations, only axial displacements and approximate knowledge of the lateral conditions at the boundaries may be sufficient to characterize the medium. In more complex cases, however, uti- lizing both axial and lateral displacements at the boundaries may be needed to successful estimate all the parameters. 108 Chapter 4. Dynamic Finite Element Inverse Problem Width (mm) D ep th (m m) Estimated Elasticity (kPa) 0 10 20 30 40 0 10 20 30 40 15 20 25 (a) Width (mm) D ep th (m m) Estimated Viscosity (Pas) 0 10 20 30 40 0 10 20 30 40 5 15 25 (b) Figure 4.8: The effect of using an inaccurate model in the inverse problem for (a) the elasticity and (b) the viscosity estimations. The displacements were computed using a model of the damping in which the damping matrix of every element is assumed to be proportional to its stiffness matrix. The inverse problem was then solved using the proposed damping structure as it results from equation (4.19). 109 Chapter 4. Dynamic Finite Element Inverse Problem Bmode (a) Width (mm) D ep th (m m) Elasticity (kPa) 10 20 30 10 20 30 40 20 40 60 80 100 (b) Width (mm) D ep th (m m) Viscosity (Pas) 10 20 30 10 20 30 40 50 100 150 200 250 (c) Width (mm) D ep th (m m) Relaxation−time (ms) 10 20 30 10 20 30 40 0.5 1 1.5 2 2.5 (d) Figure 4.9: (a) B-mode image of the gelatin phantom with a small PVA sponge inclusion. The approximate location of the inclusion is delineated in the image. The results of the (b) young’s modulus and (c) viscosity esti- mations are shown. The relaxation-time (d) image is produced as the ratio between the Young’s modulus and viscosity. The location of the inclusion and the relative values of the parameters inside and outside the inclusion are consistent with the B-mode observations and rheometry measurements. The artifacts in the images may be associated with the lack of information about the lateral boundary conditions. The images have been filtered and up-sampled for a better illustration. 110 References [1] J. Ophir, I. Cespedes, B. Garra, H. Ponnekanti, Y. Huang, and N. Mak- lad, “Elastography: ultrasonic imaging of tissue strain and elastic mod- ulus in vivo,” European Journal of Ultrasound, vol. 3, no. 1, pp. 49–70, Jan. 1996. [2] C. Pellot-Barakat, M. Sridhar, K. K. Lindfors, and M. F. Insana, “Ul- trasonic elasticity imaging as a tool for breast cancer diagnosis and research,” Current Medical Imaging Reviews, vol. 2, no. 1, pp. 157–164, Feb. 2006. [3] R. Sinkus, M. Tanter, S. Catheline, J. Lorenzen, C. Kuhl, E. Son- dermann, and M. Fink, “Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography.” Magnetic Reso- nance in Medicine, vol. 53, no. 2, pp. 372–387, Feb 2005. [4] R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink, “Viscoelastic shear properties of in vivo breast lesions measured by MR elastography,” Magnetic Resonance Imaging, vol. 23, no. 2, pp. 159–65, 2005. [5] Y. Liu, L. Sun, and G. Wang, “Tomography-based 3-d anisotropic elas- tography using boundary measurements,” IEEE Transactions on Med- ical Imaging, vol. 24, no. 10, pp. 1323 – 1333, 2005. [6] M. Sridhar, J. Liu, and M. F. Insana, “Viscoelasticity imaging using ultrasound: parameters and error analysis” Physics in Medicine and Biology, vol. 52, no. 9, pp. 2425–2443, May 2007. [7] D. D. Steele, T. L. Chenevert, A. R. Skovoroda, and S. Y. Emelianov, “Three-dimensional static displacement, stimulated echo NMR elas- ticity imaging”, Physics in Medicine and Biology, vol. 45, no. 6, pp. 1633–1648, Jun 2000. [8] H. Eskandari, S. Salcudean, and R. Rohling, “Viscoelastic parameter estimation based on spectral analysis,” IEEE Transactions on Ultrason- 111 Chapter 4. Dynamic Finite Element Inverse Problem ics, Ferroelectrics and Frequency Control, vol. 55, no. 7, pp. 1611–1625, July 2008. [9] M. Insana, J. Liu, M. Sridhar, and C. Pellot-Barakat, “Ultrasonic me- chanical relaxation imaging and the material science of breast cancer,” in Ultrasonics Symposium, 2005 IEEE, vol. 2, 18-21 Sept. 2005, pp. 739–742. [10] E. Turgay, S. E. Salcudean, and R. Rohling, “Identifying the mechan- ical properties of tissue by ultrasound strain imaging,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 221–235, 2006. [11] K. R. Raghavan and A. E. Yagle, “Forward and inverse problems in elasticity imaging of soft tissues,” IEEE Transactions on Nuclear Sci- ence, vol. 41, no. 4 pt 1, pp. 1639–1645, 1994. [12] C. Sumi, A. Suzuki, and K. Nakayama, “Estimation of shear modulus distribution in soft tissue from strain distribution,” IEEE Transactions on Biomedical Engineering, vol. 42, no. 2, pp. 193–202, Feb. 1995. [13] A. R. Skovoroda, S. Y. Emelianov, and M. O’Donnell, “Tissue elasticity reconstruction based on ultrasonic displacement and strain images,” IEEE transactions on ultrasonics, ferroelectrics, and frequency control, vol. 42, no. 4, pp. 747–765, 1995. [14] J. Bishop, A. Samani, J. Sciarretta, and D. B. Plewes, “Two- dimensional mr elastography with linear inversion reconstruction: methodology and noise analysis.” Phys Med Biol, vol. 45, no. 8, pp. 2081–2091, Aug 2000. [15] Y. Zhu, T. J. Hall, and J. Jiang, “A finite-element approach for young’s modulus reconstruction,” IEEE Transactions on Medical Imag- ing, vol. 22, no. 7, pp. 890–901, Sep. 2003. [16] E. Park and A. M. Maniatty, “Shear modulus reconstruction in dynamic elastography: time harmonic case.” Phys Med Biol, vol. 51, no. 15, pp. 3697–3721, Aug 2006. [17] F. Kallel and M. Bertrand, “Tissue elasticity reconstruction using linear perturbation method,” IEEE Transactions on Medical Imaging, vol. 15, no. 3, pp. 299–313, Jun 1996. 112 Chapter 4. Dynamic Finite Element Inverse Problem [18] D. Fu, S. F. Levinson, S. M. Gracewski, and K. J. Parker, “Non-invasive quantitative reconstruction of tissue elasticity using an iterative forward approach,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1495– 509, Aug 2000. [19] M. I. Miga, “A new approach to elastography using mutual information and finite elements,” Physics in Medicine and Biology, vol. 48, no. 4, pp. 467–80, 02/21 2003. [20] 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, no. 13, pp. 2955–2974, Jul. 2004. [21] M. M. Doyley, E. E. V. Houten, J. B. Weaver, S. Poplack, L. Duncan, F. Kennedy, and K. D. Paulsen, “Shear modulus estimation using par- allelized partial volumetric reconstruction.” IEEE Trans Med Imaging, vol. 23, no. 11, pp. 1404–1416, Nov. 2004. [22] A. S. Khalil, R. C. Chan, A. H. Chau, B. E. Bouma, and M. R. K. Mofrad, “Tissue elasticity estimation with optical coherence elastogra- phy: Toward mechanical characterization of in vivo soft tissue,” Annals of Biomedical Engineering, vol. 33, no. 11, pp. 1631–1639, 2005. [23] P. E. Barbone and A. A. Oberai, “Elastic modulus imaging: some exact solutions of the compressible elastography inverse problem.” Phys Med Biol, vol. 52, no. 6, pp. 1577–1593, Mar 2007. [24] P. E. Barbone and J. C. Bamber, “Quantitative elasticity imaging: what can and cannot be inferred from strain images.” Phys Med Biol, vol. 47, no. 12, pp. 2147–2164, Jun 2002. [25] P. E. Barbone and N. H. Gokhale, “Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions,” Inverse Problems, vol. 20, no. 1, pp. 283–96, Feb 2004. [26] J. R. McLaughlin and J.-R. Yoon, “Unique identifiability of elastic parameters from time-dependent interior displacement measurement,” Inverse Problems, vol. 20, no. 1, pp. 25–45, 2004. [27] H. Eskandari and S. E. Salcudean, “Robust estimation of tissue elas- ticity using dynamic finite elements and spectral averaging,” in 29th Canadian Medical and Biological Engineering Conference, Vancouver, Canada, June 2006. 113 Chapter 4. Dynamic Finite Element Inverse Problem [28] W. F. Walker, F. J. Fernandez, and L. A. Negron, “A method of imag- ing viscoelastic parameters with acoustic radiation force.” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1437–1447, Jun 2000. [29] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity and viscosity from measurement of shear wave speed dispersion,” Journal of the Acoustical Society of America, vol. 115, no. 6, pp. 2781–2785, 2004. [30] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelka- ram, and J. Culioli, “Measurement of viscoelastic properties of homo- geneous soft solid using transient elastography: An inverse problem ap- proach,” Journal of the Acoustical Society of America, vol. 116, no. 6, pp. 3734–3741, 2004. [31] J. Bercoff, M. Tanter, M. Muller, and M. Fink, “The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force,” IEEE transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 51, no. 11, pp. 1523–1536, 2004. [32] S. Girnyk, A. Barannik, E. Barannik, V. Tovstiak, A. Marusenko, and V. Volokhov, “The estimation of elasticity and viscosity of soft tissues in vitro using the data of remote acoustic palpation,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 211–219, 2006. [33] H. Eskandari, S. Salcudean, and R. Rohling, “Tissue strain imaging using a wavelet transform-based peak search algorithm,” IEEE Trans- actions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 54, no. 6, pp. 1118–1130, June 2007. [34] C. Gontier, M. Smail, and P. Gautier, “A time domain method for the identification of dynamic parameters of structures,” Mechanical Sys- tems and Signal Processing, vol. 7, no. 1, pp. 45–56, Jan 1993. [35] D. F. Pilkey, G. Park, and D. J. Inman, “Damping matrix identification and experimental verification,” in Proceedings of SPIE, vol. 3672, 1999, pp. 350–357. [36] S. Adhikari, “Lancaster’s method of damping identification revisited,” Journal of Vibration and Acoustics, Transactions of the ASME, vol. 124, no. 4, pp. 617–627, 2002. [37] L. E. Malvern, Introduction to the mechanics of a continuous medium, 1st ed. New Jersey: Prentice-Hall Inc., 1969. 114 Chapter 4. Dynamic Finite Element Inverse Problem [38] O. C. Zienkiewicz and R. L. Taylor, The finite element method, 5th ed. Oxford, UK: Butterworth-Heinemann, 2000. [39] R. D. Cook, D. S. Malkus, and M. E. Palesha, Concepts and Applica- tions of Finite Element Analysis, 3rd ed. New York: John Wiley & Sons, Inc., 1989. [40] M. Bro-Nielsen, “Finite element modeling in surgery simulation,” Pro- ceedings of the IEEE, vol. 86, no. 3, pp. 490–503, March 1998. [41] J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. New York: Springer, 2006. [42] D. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM Journal on Applied Mathematics, vol. 11, no. 2, pp. pp. 431–441, June 1963. [43] M. M. Doyley, P. M. Meaney, and J. C. Bamber, “Evaluation of an iterative reconstruction method for quantitative elastography,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1521–40, Jun. 2000. [44] A. A. Oberai, N. H. Gokhale, and G. R. Feijoo, “Solution of inverse problems in elasticity imaging using the adjoint method,” Phys Med Biol, vol. 19, no. 2, pp. 297–313, 2003. [45] T. Varghese and J. Ophir, “An analysis of elastographic contrast-to- noise ratio.” Ultrasound in Medicine and Biology, vol. 24, no. 6, pp. 915–924, July 1998. [46] J. F. Semblat, “Rheological interpretation of Rayleigh damping,” Jour- nal of Sound and Vibration, vol. 206, no. 5, pp. 741–744, Oct 1997. [47] R. Zahiri-Azar and S. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, Oct. 2006. [48] L. Ljung, System identification : theory for the user, 2nd ed. NJ: Prentice Hall PTR, 1999. [49] J. Ophir, S. K. Alam, B. Garra, F. Kallel, E. Konofagou, T. Korouskop, T. Varghese, “Elastography: ultrasonic estimation and imaging of the elastic properties of tissues”, Proc. Inst. Mech. Eng. [H], vol. 213, no. 3, pp. 203–233, May 2007. 115 Chapter 4. Dynamic Finite Element Inverse Problem [50] 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 and comparison with gold standard testings”, IEEE Transactions on Ul- trasonics, Ferroelectrics and Frequency Control, vol. 54, no. 3, pp. 498– 509, Feb 2007. 116 Chapter 5 Model Validation for Longitudinal Wave Propagation in a Viscoelastic Environment 4 5.1 Introduction The propagation of the wave in soft tissues is governed by local mechanical properties. In an inverse problem scheme, the unknown mechanical para- meters can be identified by analyzing the dynamic motion of the medium. While there are two types of wave that can propagate in a viscoelastic en- vironment, shear waves have received a significantly more attention than longitudinal waves in the field of tissue characterization. For a shear or transverse wave, particle displacements are in a direction perpendicular to the direction of the wave propagation, while for a longitudinal wave, the particle motions are in the same direction as the wave propagates. Shear waves can be produced inside the tissue either by external transverse motion of the boundary [1,2], or through remotely induced acoustic radiation force impulse (ARFI) [3–6]. The resulting shear wave can be tracked by ultra- sound or by magnetic resonance elastography (MRE) and the shear modulus can be reconstructed from the velocity of the wave front [7]. Longitudinal waves are assumed to have a very high propagation speed which makes them untraceable by ultrasound or low frame-rate data acquisi- tion methods. Therefore, the effect of longitudinal wave is usually eliminated by applying the curl operator to the displacement field [8]. We have recently shown that it is possible to track the compressional motion of the medium 4A version of this chapter has been submitted for a journal publication. H. Eskandari, A. Baghani, S. E. Salcudean and R. Rohling, “The Influence of the Boundary Conditions on Longitudinal Wave Propagation in a Viscoelastic Medium”. Physics in Medicine and Biology, submitted in Mar. 2009 117 Chapter 5. Longitudinal Wave Propagation using ultrasound vibro-elastography [9, 10]. Observations of wave propaga- tion in tissue-mimicking phantoms indicate that a compression travels at a much lower speed in finite media than the speed of ultrasound, which makes it possible to track the wave front by means of ultrasound. Similar results were also obtained in [11] and [12], where the relaxation behavior of phan- toms and breast tissue were observed and quantified using a compressional scheme. In [13], a method for generating low-speed longitudinal waves was pro- posed. A low-frequency compressional excitation was applied by two rods that were placed on both sides of an ultrasound probe. As a result, a wave could be seen in front of the probe, having the same speed as the shear wave. This phenomenon was explained in [13] and [14] as being the effect of superposition of the shear wave fronts of the two exciters, superimposing in the middle line and producing a longitudinally polarized shear wave through mode conversion. A Green’s function analysis in [15] revealed the effects of a coupling term in observing waves that propagate longitudinally at the speed of shear wave. In this paper we investigate the effect of boundary conditions and ex- citation on the speed of longitudinal waves in a finite medium. It will be shown by analyzing the three-dimensional (3D) wave equation and through finite element simulations that a compressional wave may travel at very high speed or at speeds as low as that of the shear wave. While a plane p-wave in an incompressible medium will travel at a speed close to that of ultrasound, depending on the boundary conditions, a slow longitudinal wave may also be generated in the medium with the main component of its particle motion in the direction of propagation. As a result, besides shear wave excitation, tracking the longitudinal waves in ultrasound elastography can be a used as another method to investigate the viscoelastic behavior of the materi- als. A 3D FEM model was adopted as a reference for the analysis and to validate the theoretical results, while the error produced by simplified 2D FEM or even 1D models are explored. Transient and harmonic simulations were performed. A linear viscoelastic medium was assumed and the effect of viscosity on the displacement spectrum is discussed. 5.2 Governing Equations in a Linear Viscoelastic Medium In the Cartesian coordinate, the three dimensional (3D) displacement vector ~u = [ux, uy, uz] T and the strain tensor ǫ = [ǫxx, ǫyy, ǫzz, γxy, γxz, γyz] T 118 Chapter 5. Longitudinal Wave Propagation have the following relationship: ǫ = ∂x 0 0 0 ∂y 0 0 0 ∂z ∂y ∂x 0 ∂z 0 ∂x 0 ∂z ∂y ~u . (5.1) In a linear viscoelastic medium, the stress tensor σ is related to the time-dependent strain tensor through a combination of the Lamé constants λ and µ, and viscosity constants λ′ and µ′. σ = λ + 2µ λ λ 0 0 0 λ λ + 2µ λ 0 0 0 λ λ λ + 2µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ ǫ + λ′ + 2µ′ λ′ λ′ 0 0 0 λ′ λ′ + 2µ′ λ′ 0 0 0 λ′ λ′ λ′ + 2µ′ 0 0 0 0 0 0 µ′ 0 0 0 0 0 0 µ′ 0 0 0 0 0 0 µ′ ∂tǫ , (5.2) where ∂t is the time-derivative. Here, µ and µ ′ are often referred to as shear elasticity and shear viscosity, and λ′ is known as the elongational or compressional viscosity. The equation of motion for a body that is subjected to an external force field ~f , is: ρ ∂2t ~u = ∂xσxx + ∂yσxy + ∂zσxz∂xσxy + ∂yσyy + ∂zσyz ∂xσxz + ∂yσyz + ∂zσzz + ~f , (5.3) where ~u is the displacement vector as a function of time and position, and ρ is the density. By inserting from (5.1) and (5.2) into (5.3), a hyperbolic partial differential equation (PDE) results which describes the wave propagation in 119 Chapter 5. Longitudinal Wave Propagation a linear viscoelastic medium [18]: ρ ∂2t ~u = (λ+ µ)∇∇ · ~u+ µ∇2~u+ (λ′ + µ′)∇∇ · (∂t~u) + µ′∇2(∂t~u) + ~f . (5.4) The gradient and divergence operators are denoted by (∇) and (∇·), respec- tively. The vector Laplacian ∇2 is defined as the divergence of the gradient of the vector. Note that in the Cartesian coordinates ∇· ~u = ǫxx+ ǫyy+ ǫzz, which is a measure of the volumetric change in the medium and is called the dilatation. The Young’s modulus E is defined as: E = (1− 2ν)λ+ 2µ , (5.5) where ν is the Poisson’s ratio, which takes values between 0 and 0.5. A value close to 0.5 denotes incompressibility of the material as is the case for soft tissues. Different assumptions have been made in the literature to re- duce the number of viscosity parameters. For example, if the viscous effects are assumed to follow the elastic effects in the material, then λ′/λ = µ′/µ. However, there has been little justification for such similar behavior of visc- sity and elasticity which makes this assumption implausible [16]. Another assumption as made by [4] and [8] is that the compressional viscosity λ′ is negligible compared to the shear viscosity. To our knowledge, there is no experimental validation for this assumption for incompressible soft tissue. Another approach is to translate the idea of elastic incompressibility into the viscosity domain by assuming that the bulk viscosity is negligible, so that the fluid incompressibility condition holds. When the divergence of flow is zero, the fluid is considered to be incompressible and the term asso- ciated with the bulk viscosity in the Navier-Stokes equations will vanish. In this case, the power dissipation in the medium is only diviatoric dissipation associated with the distortion or shape-change, rather than dilatational dis- sipation that involves a change of volume [17]. The bulk viscosity is defined as κ = λ′ + 2µ′/3 which can be set to zero to obtain: λ′ = −2 3 µ′ . (5.6) Bulk viscosity is known to be negligible compared to shear viscosity for rubber-like materials as well as some liquids and gases [18]. We proceed with the assumption that both conditions of elastic incompressibility (ν ≈ 0.5 in (5.5)) and fluid incompressibility (5.6) are marginally satisfied in soft tissue. 120 Chapter 5. Longitudinal Wave Propagation 5.3 Wave Equation and the Propagation Speed In a purely elastic medium, the wave equation in (5.4) becomes: ρ ∂2t ~u = (λ+ µ)∇∇ · ~u+ µ∇2~u+ ~f . (5.7) The general solution of equation (5.7) can be decomposed into a divergence- free and a curl-free component [18]. The divergence-free part tends to main- tain the volume and is referred to as the shear or s-wave which travels at a speed of: cs = √ µ ρ . (5.8) 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 = √ λ+ 2µ ρ . (5.9) While pure shear or dilatational waves travel at the designated speeds mentioned above, the influence from the boundary conditions may give rise to waves that travel at other velocities. A few examples are presented in the following sections. 5.3.1 Longitudinal Wave in a Finite Medium Approximating a Thin Rod Let us first consider the case where the excitation is mainly compressional in the x direction. If the exciter is large relative to the cross-section such that a plane-wave is induced in the medium, and the sides of the region have free boundaries, all the stress components except for σxx, may be ignored. Note that for a viscoelastic case, the strain distribution will be affected by the viscous stresses as well as the elastic stresses, however due to the compressional excitation and the plane-wave assumption, σxx is dominant and the approximation should be valid. Using this condition in (5.2),[ (λ+ µ) + (λ′ + µ′)∂t ] σxx = ( µ+ µ′∂t ) [ (3λ+ 2µ) + (3λ′ + 2µ)∂t ] ǫxx . (5.10) For the interior region that is not experiencing any internal force, (5.3) is reduced to ∂xσxx = ρ ∂ 2 t ux. This can be substituted into (5.10) to obtain 121 Chapter 5. Longitudinal Wave Propagation the viscoelastic wave equation in the x -direction: ρ [ (λ+ µ) + (λ′ + µ′)∂t ] ∂2t ux = ( µ+ µ′∂t ) [ (3λ+ 2µ) + (3λ′ + 2µ)∂t ] ∂2xux . (5.11) This is a third order PDE with respect to the variable t. The case of longitudinal wave propagation in a viscoelastic thin rod has been addressed more than seven decades ago by Thompson [16]. The conditions described in this section in which non-axial stress components become zero approximate a thin rod situation. Therefore, following [16], by taking the Laplace transform in the x domain (L{f(x)} = F (s)) and the Fourier transform in the t domain (F{f(t)} = f̂(jω)), the characteristic equation of (5.11) becomes: −ω2ρ+ s2 (µ+ jωµ ′) [(3λ+ 2µ) + jω(3λ′ + 2µ′)] [(λ+ µ) + jω(λ′ + µ′)] = 0 . (5.12) There are three roots associated with this characteristic equation with re- spect to ω. However, considering that at relatively low frequencies (λ+µ)≫ (λ′+µ′)ω for incompressible materials, the root that corresponds to the de- nominator is very large and may be ignored [16]. Also, with the assumption of zero bulk viscosity, 3λ′ + 2µ′ = 0. Therefore, equation (5.12) reduces to: −ω2ρ+ s2 (µ+ jωµ ′)(3λ + 2µ) (λ+ µ) = 0 . (5.13) Given the incompressibility of the material and the definition of Young’s modulus, this can be written as −ω2ρ + s2(E + jωη) = 0, which results in the following PDE: ρ∂2t ux = E∂ 2 xux + η∂ 2 x∂tux . (5.14) η is called the coefficient of normal viscosity, where, η ≈ 3µ′ . (5.15) Equation (5.14) describes the propagation of the wave in a 1D environ- ment. This PDE can be discretized using the well-known distributed system of mass-spring-damper (MSD) elements in series [19]. Based on (5.15), in order for the 1D model to follow the behavior of the 3D model in its best sense, the equivalent viscosity of the 1D model should be three times the shear viscosity. This will be taken into account in the following simulations where behaviors of different models are compared. Also, given that the Young’s modulus in incompressible materials is approximately three times 122 Chapter 5. Longitudinal Wave Propagation the shear modulus, the ratio between the observed viscous effect and the elastic or shear modulus, entitled relaxation-time, would be independent of the modality of the excitation, whether it is shear or longitudinal. Equation (5.14) indicates that a longitudinal wave is present in the medium, propagating at a much lower speed than cp. In the purely elastic case, without the influence of viscosity, the speed is: c0 = √ E ρ . (5.16) This value, being approximately 1.7 times the shear wave speed for in- compressible materials, denotes that a low-speed longitudinal wave can be induced in the medium by applying certain boundary conditions. Note that if the viscosity is also taken into account, the effective modulus from (5.14) will be √ E2 + ω2η2 instead of E, which results in a speed equal to c ≈ c0(1 + ω2τ2/4), where τ = η/E is the relaxation-time. However, with the small viscosity and relaxation-time in soft tissue, this correction term is negligible at low frequencies and equation (5.16) can be used to calculate the wave speed. 5.3.2 Longitudinal Wave Induced by a Point Source Exciter Suppose a point source excitation is applied axially to a block of homo- geneous material in which radial symmetry can be assumed. In this case, it is more convenient to analyze the problem in cylindrical coordinates as shown in Fig. 5.1. The stress-strain relationship in the cylindrical coordi- nate system is the same as (5.2). Similar to the analysis in Section 5.2, in the cylindrical coordinates where ǫ = [ǫxx, ǫrr, ǫθθ, γxr, γxθ, γrθ] T and ~u = [ux, ur, uθ] T , the strain tensor is defined as follows: ǫ = ∂x 0 0 0 ∂r 0 0 1 r 1 r ∂θ ∂r ∂x 0 1 r ∂θ 0 ∂x 0 1 r ∂θ (∂r − 1r ) ~u . (5.17) The equation of motion in the x -direction in this case becomes: ρ ∂2t ux = ∂xσxx + ∂rσxr + 1 r σxr + 1 r ∂θσxθ + fx . (5.18) 123 Chapter 5. Longitudinal Wave Propagation Note that if all the shear stress components are zero, the resulting wave equation would be that of the p-wave in a viscoelastic medium. For the interior, where fx = 0, (5.18) can be expanded as: ρ ∂2t ux = [(λ+ 2µ) + (λ ′ + 2µ′)∂t] ∂x(ǫxx + ǫrr + ǫθθ) −2(µ+ µ′∂t)∂x(ǫrr + ǫθθ) +(µ+ µ′∂t)∂rǫxr + 1 r (µ+ µ′∂t)ǫxr + 1 r (µ+ µ′∂t)∂rǫxθ . (5.19) It can be shown that equation (5.19) can be simplified to: ρ ∂2t ux = [(λ+ 2µ) + (λ ′ + 2µ′)∂t] ∂x(ǫxx + ǫrr + ǫθθ) −2(µ+ µ′∂t)1r [∂r(rω̄θ)− ∂θ(ω̄r)] , (5.20) where ω̄r and ω̄θ are rotations about the r̂ and θ̂ axes: 2ω̄θ = ∂zur − ∂ruz , and, 2ω̄r = 1 r ∂θuz − ∂zuθ . (5.21) To obtain a relationship between the radial and angular strains on the axis of excitation, l’Hôpital’s rule can be utilized at r → 0. According to this theorem, if two functions are differentiable in a neighborhood of a given point while both of them converge to zero, the indeterminate limit of their quotients would be equal to the quotient of their derivatives at that point. By applying this rule to the angular strain, one gets: ǫθθ|r=0 = lim r→0 1 r ur = lim r→0 ∂rur = ǫrr|r=0 . (5.22) Also, due to radial symmetry in this case, σθθ = 0, which according to (5.2), leads to:[ (λ+ 2µ) + (λ′ + 2µ′)∂t ] ǫθθ + (λ+ λ ′∂t)(ǫxx + ǫrr) = 0 . (5.23) Inserting from (5.22) into (5.23), at the axis of excitation: ǫrr = ǫθθ = −(λ+ λ′∂t) 2 [(λ+ µ) + (λ′ + µ′)∂t] ǫxx . (5.24) Inserting (5.24) into (5.20) and also noting that all rotations would be 124 Chapter 5. Longitudinal Wave Propagation u r u è u x r è x fx Figure 5.1: Displacement vector in the cylindrical coordinate system. 125 Chapter 5. Longitudinal Wave Propagation zero on the excitation axis, the following results: ρ [ (λ+ µ) + (λ′ + µ′)∂t ] ∂2t ux = ( µ+ µ′∂t ) [ (λ+ 2µ) + (λ′ + 2µ)∂t ] ∂2xux . (5.25) With the same approach as in section 5.3.1 and given the large value of λ compared to µ, µ′ and λ′ in incompressible materials, the wave equation simplifies to: ρ∂2t ux = µ∂ 2 xux + µ ′∂2x∂tux . (5.26) Therefore, given radial symmetry in the medium, a point-source exciter will generate a longitudinal wave that travels at the shear speed: cps = √ µ ρ , (5.27) Although ω̂θ is not identically equal to zero on the axis of excitation, the assumption of relatively small ω̂θ is validated in Section 5.5 by observing the same wave as predicted above. Also, note that using (5.23), the axial stress will be: σxx = [(λ+ 2µ) + (λ ′ + 2µ′)∂t] ǫxx + (λ+ λ ′∂t)(ǫrr + ǫθθ) = (µ+µ ′∂t)[(3λ+2µ)+(3λ′+2µ′)∂t] [(λ+µ)+(λ′+µ′)∂t] ǫxx ≈ E ǫxx + η ∂tǫxx , (5.28) where η is as in equation (5.15). In other words, the apparent elastic con- stant in the axial direction is the Young’s modulus rather than the shear modulus, however the velocity of the wave is determined by the shear mod- ulus according to (5.27). 5.4 Determining the Wave Speed from the Spectra of the Displacements In this section, the effect of the viscosity on the wave equation and in partic- ular on the displacement spectrum will be analyzed. It will be shown that the wave speed can be estimated from the spectrum of the displacements in a viscoelastic medium. The result from this section will be used to estimate the wave speed from the simulation results in Section 5.5, where the effect of the boundary conditions on the wave speed is explored. As shown in [10], if the medium is bounded from one end, the solution to (5.14) in the frequency 126 Chapter 5. Longitudinal Wave Propagation domain is in the following form: û(x, ω) = A0 [cos(βx)sinh(αx) + jsin(βx)cosh(αx)] , (5.29) where û is the Fourier transform of ux, A0 is the amplitude of the vibration at frequency ω and α and β can be calculated by taking the Laplace transform of (5.14) in the x domain and the Fourier transform in the t domain [10]: (E + jωη)s2 + ω2ρ = 0 . (5.30) The relaxation-time is defined as τ = η/E. Given that in soft tissue, ωτ ≪ 1 for ω < 100Hz, the solution to (5.30) can be written as: s = ±(α+ jβ) = ± jω c0 √ 1 + jωτ ≈ ±jω c0 (1− jωτ 2 ) , (5.31) where c0 is defined in (5.16). Hence, α ≈ ω 2τ 2c0 , β ≈ ω c0 . (5.32) 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 (5.29) is: |û(x, ω)| = sinh 2(αx) + sin2(βx) sinh2(αL) + sin2(βL) . (5.33) Equation (5.33) is a 2D function of ω and x, from which the wave speed c0 can be obtained. Two points on this functions are of particular interest to us: (i) The frequency f1 at which the displacement at x = L/2 is larger than other locations. (ii) The frequency f2 at which the displacement at x = L/2 has its maxi- mum value compared to other frequencies. Note that for a purely elastic medium, these two points will be identical. So that, the wave speed can be estimated by finding either the frequency that maximizes the displacement at the middle point or the frequency for 127 Chapter 5. Longitudinal Wave Propagation which the middle point has the highest displacement magnitude compared to other locations. Figure 5.2 shows the displacement magnitude in equation (5.33) as a function of frequency and depth for a low viscosity and a high viscosity medium. A Young’s modulus of 10kPa and density of 1000kg/m3 were assumed, thus c0 = 3.16m/3 and f0 = c0/(2L) = 26.4Hz. First, with a viscosity of 5Pas, α and β were calculated from (5.31). The displacement magnitude and the equal contour plots are shown in Figs. 5.2(a) and 5.2(b) where the frequencies as noted by case (i) and case (ii) are shown as f1 and f2 respectively. It can be seen that if the viscosity is very low, the two cases are almost the same and both f1 and f2 occur very close to f0. In the second case, a viscosity of 30Pas was assumed and the displacement magnitudes were calculated. The magitude and contour plots are depicted in Figs. 5.2(c) and 5.2(d). In this case, f1 = 33.2Hz and f2 = 24.4Hz, therefore, using f2 to reconstruct the wave speed would be more accurate. In the following, we will explore how the wave speed can be reconstructed from case (i) using f1 and from case (ii) using f2. The errors that each approximation produces will be quantified. The better method will be adopted for wave speed analysis in the remaining sections of the paper. 5.4.1 Case (i) According to (i), ω1 should be found such that, ∂x |û(x, ω1)| = 0, at x = L/2. By differentiating (5.33) with respect to x, ∂x |û(x, ω)| = αsinh(2αx) + βsin(2βx) sinh2(αL) + sin2(βL) . (5.34) This equation should be zero at [x, ω] = [L/2, ω1]. Therefore, 1 2 ω1τsinh(αL) + sin(βL) = 0 . (5.35) Since we assumed that ωτ ≪ 1 and L≪ 1, therefore, sinh(αL) ≈ αL. As a consequence, it can be shown from (5.35) that sin(βL) ≈ 0 while βL ≫ 0. Thus, sin(βL) ≈ π − βL. From (5.35) and (5.32): ω31τ 2L 4c0 + π − ω1L c0 = 0 . (5.36) If ω1 = 2πf1, the wave speed can be estimated as follows: c0 = 2f1L ( 1− π2f21 τ2 ) . (5.37) 128 Chapter 5. Longitudinal Wave Propagation Frequency (Hz) D ep th (c m) 10 20 30 40 50 60 0 1 2 3 4 5 6 10 20 30 40 50 (a) Frequency (Hz) D ep th (c m) 20 40 60 6 5 4 3 2 1 0 f1, f2 (b) Frequency (Hz) D ep th (c m) 10 20 30 40 50 60 0 1 2 3 4 5 6 0.5 1 1.5 (c) Frequency (Hz) D ep th (c m) 20 40 60 6 5 4 3 2 1 0 f1f2 (d) Figure 5.2: The magnitude of the displacement in a 1D model of the wave equation with the purely elastic wave speed of 3.16m/s, length of 6cm, and Young’s modulus of 10kPa. (a) and (b) are the magnitude and contour plots for the case that 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 case that the viscosity of the medium is 30Pas. In this case, f1 over-estimates the actual value, while f2 results in a relatively small under-estimation. 129 Chapter 5. Longitudinal Wave Propagation Apparently in this case, if one uses 2f1L instead of (5.37), the wave speed will be over-estimated by a small fraction. 5.4.2 Case (ii) In this case, ω2 should be found such that ∂ω |û(x, ω)| = 0 at [x, ω] = [L/2, ω2]. By taking the derivative of (5.33) with respect to ω and putting it equal to zero for x = L/2, the following will be obtained: (ω2τ sinh(αL) + sin(βL)) ( sinh2(αL) + sin2(βL) ) = 2 (ω2τ sinh(2αL) + sin(2βL)) ( sinh2( αL 2 ) + sin2( βL 2 ) ) . (5.38) Using the characteristics of the trigonometric and hyperbolic functions, it can be shown that: sinh2(αL) + sin2(βL) = (cosh(αL) + cos(βL)) (cosh(αL) − cos(βL)) , (5.39) and, sinh2( αL 2 ) + sin2( βL 2 ) = cosh(αL)− cos(βL) 2 . (5.40) Inserting (5.39) and (5.40) into (5.38) gives: (ω2τ sinh(αL) + sin(βL)) (cosh(αL) + cos(βL)) = ω2τ sinh(2αL) + sin(2βL) , (5.41) which can be simplified to: sin(βL) = ω2τ sinh(αL) . (5.42) Using the same assumptions as noted in Case (i), the wave speed can be calculated as: c0 = 2f2L ( 1 + π2f22 τ 2 ) . (5.43) where ω2 = 2πf2. If a small amount of under-estimation is tolerable, the wave speed can be approximated as 2f2L. By comparing (5.37) and (5.43), it can be seen that f2 < f1, therefore the second-order term in (5.43) is smaller than the second-order term in (5.37). As a result, it would be a more accurate estimation to approximate the wave speed as 2f2L rather than 2f1L, should it be necessary, for example if the relaxation-time of the medium is unknown. Also, note that although 130 Chapter 5. Longitudinal Wave Propagation equation (5.14) has been used here to analyze the accuracy of the two cases above, the same argument is valid for other forms of the wave equation, for example when shear modulus replaces Young’s modulus as in (5.26). 5.5 Results 5.5.1 Dynamic Finite Element Analysis Numerical analysis of the three-dimensional (3D) wave equation is possible through discretization of equation (5.4) using Ritz-Galerkin method [20]. The resulting finite element model can be used to explain the static and dy- namic deformations in the medium in a discretized platform. An equation can be obtained to describe the relationship between the nodal displace- ments, boundary conditions and local material properties: Ku(t) +Bu̇(t) +Mü(t) = f(t) , (5.44) where K, B and M are the stiffness, viscosity and mass matrices of the elements, u(t) is the nodal displacement vector and f(t) is a vector contain- ing the boundary conditions. Equation (5.44) should hold at all times. This equation can be written in the frequency domain as follows: (K + jωB − ω2M)û = f̂ . (5.45) A method to generate an accurate viscosity matrix based on a linear viscoelastic model of the continuum has been proposed in [21]. An implicit solution for the transient case of (5.44) with step-size ∆t can be obtained as [22]: K̄u(t+∆t) = f̄(t) , K̄ =K + B ∆t + M ∆t2 , f̄(t) = f(t) + 2M ∆t2 u(t) + ( B 2∆t − M ∆t2 )u(t−∆t) . (5.46) A finite element simulation can be used to verify the theory presented in Section 5.3. In particular, numerical tests can facilitate analysis of the propagation of the fast and slow longitudinal waves in soft materials. For this purpose, a region of size 6×4×4(cm3) was meshed by eight-node brick elements, having a total of 25 × 7 × 7 nodes. The bottom surface at x = 0 was confined from axial movement in the x direction and the excitation was 131 Chapter 5. Longitudinal Wave Propagation X Z Y 4cm 6 cm 4 cm[0, 0, 0]T Excitation Figure 5.3: The geometry of the 3D region used for FEM simulations. 132 Chapter 5. Longitudinal Wave Propagation applied axially at the top surface, as shown in Fig. 5.3. Although some of the derivations in section 5.3 have been based on radial symmetry or thin rod assumptions, the Cartesian-based FEM model utilized here does not satisfy those conditions. However, the geometry used here can additionally provide an insight into the longitudinal wave propagation in a non-ideal situation. In this manuscript, the middle line refers to the line in the phantom with {y, z} = [2, 2]cm, the middle point or middle node refers to the point or node on top of the phantom with {x, y, z} = [6, 2, 2]cm, and the center point is the point at {x, y, z} = [3, 2, 2]cm. 5.5.2 Simulation of the Wave Propagation The 3D region in Fig. 5.3 was assumed to consist of a homogeneous material with a Young’s modulus of 10kPa, a viscosity of 10Pas and a density of 1000kg/m3 . A Poisson’s ratio of ν = 0.495 is often considered for finite element analysis of soft tissues [23]. With the given Young’s modulus and Poisson’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 Section 5.3 for the following cases, when: (a) The excitation was applied to all of the nodes at the top surface of the phantom. The four vertical sides of the phantom were free, thus bulging was allowed. (b) The excitation was applied to all of the nodes at the top surface of the phantom. The four vertical sides of the phantom were confined with slip boundaries such that bulging was not allowed. (c) The excitation was applied to a single node at the center of the top surface, i.e. the middle node. The four vertical sides of the phantom were free, thus bulging was allowed. In all these three cases, a slip boundary condition was imposed on the bottom of the phantom, while being axially confined. In case (a), an elastic longitudinal wave is expected to be generated in the medium, traveling at a speed designated by (5.16). In case (b), a volume change is not possible due to the boundary conditions. Therefore, a dilatational wave will be generated at a speed determined by (5.9). Finally, in case (c), a point source compressional exciter is modeled within the resolution of the finite element mesh. As shown in section 5.3.2, the speed of longitudinal component of the 133 Chapter 5. Longitudinal Wave Propagation generated wave can be approximated by (5.27). These theoretical derivations will be tested in the following simulations. For each case, a frequency sweep analysis has been performed through solving equation (5.45) at a range of frequencies. The transfer function between the displacements at the center point and the middle point was then measured by dividing the displacements spectra at the two locations. As explained in section 5.4, by determining the frequency at which the peak of the magnitude of this transfer function occurs, the wave speed can be estimated given the small relaxation-time of the medium. The transfer function magnitudes for the three cases above are shown in Fig. 5.4(a) and the phases are plotted in Fig. 5.4(b). The peak of the three transfer functions occur at 24.4Hz, 153Hz and 14.6Hz for cases (a), (b) and (c), respectively. These frequencies can be multiplied by a wavelength of 12cm, which is twice the axial extent of the phantom, to obtain the longitudinal wave speed. Therefore, the estimated wave speed for these cases would be respectively, 2.93m/s, 18.36m/s and 1.75m/s. The actual values as computed from equations (5.9), (5.16) and (5.27) are 3.16m/s, 18.38m/s and 1.82m/s. Under different loading conditions, the effective viscosity of the medium would change. In case (a), the coefficient of viscosity can be calculated from (5.15) to be 30kPa, which should be divided by the Young’s modulus to get a relaxation-time of τa = 1ms. In case (b), the effective viscous modulus would be λ′+2µ′, which is a counterpart of effective elastic modulus λ+2µ. Given a bulk viscosity of zero, λ′ + 2µ′ = 4µ′/3 = 13.3Pas. Thus, the relaxation-time in this case would be τb = 39µs. Note that in reality where λ is in the order of 109, the viscosity for p-wave propagation would be even less significant. Lastly, for case (c), the viscosity would be equal to µ′ according to (5.27), which can be divided by the shear modulus to get a relaxation-time of τc = 1ms. Having the relaxation-time of the medium, the estimation of the wave speed can be enhanced by using equation (5.43). The resulting approximations of the longitudinal wave speed are, 2.95m/s, 18.37m/s and 1.76m/s for the three cases, respectively. The presence of viscosity has a crucial effect on the peak value of the transfer function magnitude. Higher viscosity results in higher dissipation and less power transfer in the medium. Figure 5.5 shows the effect of chang- ing the viscosity on the magnitude of the transfer function for the same model as explained above. The elastic wave was generated in this case by compressing the top of the phantom and allowing the sides to move freely as in case (a). Note that the location of the peak and thus the speed of propagation do not change significantly, while the attenuation depends on 134 Chapter 5. Longitudinal Wave Propagation 101 102 10−3 10−2 10−1 100 101 Frequency (Hz) |H( f)| TF Magnitude at the Center Point Free sides Bounded sides Point source (a) 101 102 −6 −4 −2 0 Frequency (Hz) ∠ H (f) [ra d] TF Phase at the Center Point Free sides Bounded sides Point source (b) Figure 5.4: The magnitude (a) and the phase (b) of the transfer function between the center point of the phantom and the middle point at the top. Solid line shows the case where the four vertical sides of the phantom are free to move and the dashed lines are for the case where the sides of the phantom are confined from bulging. The gray line shows the transfer function for a one-node excitation at the middle point. 135 Chapter 5. Longitudinal Wave Propagation the damping value. To show the transient behavior of the model, a step displacement with 1mm amplitude was applied to the top of the simulated phantom. Figure 5.6 shows how the step propagates along the middle line of the phantom. The previous three cases were explored here as well. First, an elastic wave was generated in the unbounded medium by applying the excitation to the top surface. Based on the peak of the displacement waveform in Fig. 5.6(a), it takes approximately 21ms for the step to travel 6cm and reach the bottom. The speed is therefore estimated as 2.9m/s. Also from Fig. 5.6(a), the resonance of the phantom in response to the step is at 24.4Hz which complies with the previous result. The magnified version of this figure is shown in Fig. 5.6(b) for a few nodes on the middle line. This figure clearly illustrates the way that information travels in the medium. It can be seen that a fast traveling wave propagates in the phantom, which reaches the bottom after 2.8ms. This corresponds to a speed of 21.4m/s which is comparable to the p-wave speed in the medium, i.e. 18.38m/s. It should be noted that determining the temporal location of the peak only gives an estimate of the wave-front rather than rendering it accurately. For the second test as in case (b), a compressional step excitation was applied to the top surface while the material was bounded from the sides, as shown in Fig. 5.6(c) and 5.6(d). Since the p-wave is the only wave that is produced in the medium, a magnified version in Fig. 5.6(d) only shows the presence of only one wavefront at a speed of 18.2m/s. Finally, the middle node on the top of the phantom was excited by a step function. The middle line displacements are depicted in Fig. 5.6(e). While the elastic wave speed is estimated at 2.1m/s, a p-wave can be identified in the magnified version in Fig. 5.6(f) with a speed of 21.4m/s. Note that in all cases, the resonances that are produced in the phantom are coherent with the peaks of the transfer functions in Fig 5.4(a). Table 5.1 gives the theoretical as well as the estimated wave speeds. The estimated values that were calculated based on the peak of the transfer function for the center node, before and after viscosity compensation are shown. Also, the wave speed that is recovered from the temporal peak of the wave-front are reported for the three cases of boundary conditions. A sample transient wave propagation for different sizes of exciter is shown in Fig. 5.7. The boundary conditions were similar to that in case (a) and (c), while the exciter size varies. In the left column, the exciter is the smallest possible, modeling a point-source within the resolution of the FEM mesh. The exciter dimensions were increased linearly from left to right until it covered the top surface of the phantom. The snapshots are taken at 136 Chapter 5. Longitudinal Wave Propagation 101 102 10−1 100 101 Frequency (Hz) |H( f)| Effect of Viscosity on the Magnitude of TF 1 Pas 5 Pas 10 Pas 20 Pas 100 Pas Figure 5.5: Effect of the viscosity on the magnitude of the transfer function between the center point and the middle point. The excitation was applied to the entire top surface while the sides were free to bulge. Theory Estimation (m/s) (m/s) Spectrum Spectrum with Transient Excitation peak visc. compen. step (a) Unbounded 3.16 2.93 2.95 2.9 (b) Bounded 18.38 18.36 18.37 18.2 (c) Point source 1.83 1.75 1.76 2.1 Table 5.1: The theoretical and estimated wave speed for three different conditions. The wave speed, estimated with three methods are shown: first, from the frequency at which the spectrum had a peak; second, the previous value was enhanced to account for the effect of viscosity; third, from the transient wave-front due to a step excitation 137 Chapter 5. Longitudinal Wave Propagation 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 Time (ms) D is pl ac em en t (m m) Unbounded Elastic Excitation (a) 0 5 10 15 20 25 Time (ms) D is pl ac em en t L og ar ith m Unbounded Elastic Excitation Excitation Second node Center node Last node Bottom node (b) 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 Time (ms) D is pl ac em en t (m m) Bounded Medium (c) 0 5 10 15 20 25 0 0.2 0.4 0.6 0.8 1 Time (ms) D is pl ac em en t (m m) Bounded Medium (d) 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 Time (ms) D is pl ac em en t (m m) Point−source Excitation (e) 0 5 10 15 20 25 Time (ms) Point−source Excitation D is pl ac em en t L og ao rit hm (f) Figure 5.6: Displacements at the middle line for a step excitation. The left column has a larger time-span to show the low-speed wave, where each line represents the displacement of one node. The high-speed wave can be seen in the short time-span of the right column. In (a) and (b) the excitation was applied to the top surface when the four sides were free. (c) and (d) show the response for a single-node exciter at the middle node. In (e) and (f) the medium was bounded and could not bulge. Note that the displacement axis in (b) and (d) are in a logarithmic scale. 138 Chapter 5. Longitudinal Wave Propagation different times from the middle plane, where z = 2cm from Fig. 5.3. The dependency of the wave profile, speed and attenuation on the size of the exciter is noticeable, especially for the left-most column. 5.5.3 Model Comparison Generally, model simplification produces systematic errors in the analysis. The error of using a certain simplified model may drastically change the behavior of the system in some ways, but not imposing a significant change in other respects. For soft tissues, it is not possible to derive an ideal model that predicts all the linear and nonlinear behaviors within the material. For small deformations, a linear 3D dynamic FEM model is often used to describe the viscoelastic changes in the material. However, due to extensive computation requirement or lack of 3D data, the 3D FEM model can be substituted by simpler models with fewer dimensions. In this section, four models are compared: 3D FEM, 2D plane-stress, 2D plane-strain and 1D mass-spring-damper models. The geometry and the boundary conditions were the same as in the previous simulations. For the MSD and FEM models the resolution in the x direction was 25 nodes, the resolution of the FEM models in the y direction was 7 nodes and there was also 7 nodes along the z axis for the 3D model. The medium was homogeneous and the material properties were the same as those in the previous simulation. For FEM tests, the model was free at the sides and the same excitation described before was applied to the entire top part of the phantom. In the first test, the transfer function between the center point and middle point was evaluated for all of the four models. Note that in the 1D case, the entire model is supposed to represent the middle line and the transfer function is between the center node and the excited node. The magnitude and phase of the transfer functions are depicted in Fig. 5.8. In another test, the middle section of the phantom was assumed to have different viscosity and elasticity properties than the background. Counting from the bottom, elements 11 through 16 had a Young’s modulus of 30kPa and a viscosity of 70Pas. In the 2D and 3D models, this contrast was only along the x axis. For the 1D model, the viscosity associated with the dampers were assumed to be three times that of the finite element models based on (5.15). The low-frequency characteristics of these four models are shown in Fig. 5.9. The models were simulated at 5Hz and the steady-state displacements and strains were calculated. Figure 5.9 shows the magnitude and phase of the displacements as well as the strains on the middle-line of the phantom. The reciprocal of the strain magnitude can be used to 139 Chapter 5. Longitudinal Wave Propagation 1 node 6 3 0 3×3 nodes 5×5 nodes 7×7 nodes 0.2 ms 6 3 0 8.4 ms 16.8 ms 25.0 ms 6 3 0 6 3 0 6 3 0 33.4 ms 6 3 0 41.6 ms 50.0 ms Width (cm) 0 2 4 6 3 0 0 2 4 0 2 4 D ep th (c m) 0 2 4 Figure 5.7: Transient wave fronts for four different exciter sizes. The exciter has the smallest possible size at the left column, modeling a point-source. The exciter size was increased from left to right where it finally covered the top surface of the phantom. 140 Chapter 5. Longitudinal Wave Propagation 101 102 10−2 10−1 100 Frequency (Hz) |H( f)| TF Magnitude for Different Models 3D FEM Plane−strain Plane−stress 1D MSD (a) 101 102 −8 −7 −6 −5 −4 −3 −2 −1 0 1 Frequency (Hz) ∠ H (f) TF Phase for Different Models (b) Figure 5.8: Transfer functions of different models for a homogeneous block with free sides under elastic deformation. 3D FEM (solid black line), 2D plane-strain (dashed orange line), 2D plane-stress (solid red line with circle markers) and 1D mass-spring-damper (dash-dot blue line) models are com- pared. (a) is the magnitude and (b) is the phase of the transfer function between the center point and middle point of the phantom. 141 Chapter 5. Longitudinal Wave Propagation 5 10 15 20 25 0 0.2 0.4 0.6 0.8 1 Node Number D is pl ac em en t (m m) Displacement Magnitude 3D FEM Plane−strain Plane−stress 1D MSD (a) 10 12 14 16 18 0.4 0.5 0.6 0.7 Node Number D is pl ac em en t (m m) Displacement Magnitude (b) 5 10 15 20 25 −0.02 −0.01 0 0.01 0.02 0.03 Node Number Ph as e (ra d) Displacement Phase (c) 5 10 15 20 0.005 0.01 0.015 0.02 0.025 Element Number M ag ni tu de Axial Strain Magnitude (d) 5 10 15 20 −0.12 −0.08 −0.04 0 0.04 Element Number Ph as e (ra d) Axial Strain Phase (e) Figure 5.9: Comparison between four different models at 5Hz in a three- layered phantom: 3D FEM, 2D plane-strain, 2D plane-stress and 1D mass- spring-damper models. (a) shows the magnitude of the displacements at the middle line which is magnified in (b) for a closer view. The phases are shown in (c). (d) is the strain magnitude and (e) is the phase of the strain divided by the angular frequency. 142 Chapter 5. Longitudinal Wave Propagation obtain the relative elasticity of the medium and the strain phase can be divided by the angular frequency to obtain the relaxation-time [9]. The main shortcoming with the 1D model is that it does not take into account the boundary conditions, hence the estimated parameters will also contain the artifacts present in the strain images. Soft tissue deformation is often analyzed in the so-called elastography methods by using a 1D model of the wave equation or 2D FEM. In order to study the error produced by using such simplified models, a tissue-mimicking phantom was constructed using bovine skin gelatin. Figure 5.10 shows the phase of the displacements at 20Hz for a 6.5×5.5×3.5cm3 tissue-mimicking gelatin phantom with a viscous inclusion. The phantom, being 6.5cm in the x direction was imaged with ultrasound up to a depth of 4cm. The boundary conditions were the same as the simulation depicted in Fig. 5.3. To capture data from within the phantom, a Sonix RP ultrasound machine (Ultrasonix Medical Corp., Richmond, BC, Canada) was used. RF signals were collected during excitation at 1300Hz and axial displacements were estimated using a time-domain cross-correlation approach with prior estimates [24]. The control parameters and accuracy of the displacement estimates have been discussed in [9]. The excitation was a Gaussian white noise bandpass filtered at 3−30Hz. The transfer function for each block at discrete frequency values were obtained using Welch windowing algorithm [25]. The average coher- ence value at 20Hz throughout the line of acquisition was 99.99%. This frequency was chosen arbitrarily at the range of excitation, exhibiting a high coherence value. The phantom was constructed with 18% by weight of bovine skin gelatin (type G9382-B, Sigma-Aldrich Inc., Oakville, ON, Canada) in water, having a PVA sponge (Ceiba Technologies, Chandler, AZ, USA) in the middle as the inclusion. The x coordinate range of the inclusion was from 1.1cm away from the fixed ultrasound probe to 2.0cm away from it. Rheometric measurement of elasticity was not performed on the materials and an approximate value of 30kPa was assumed because of their close palpable stiffnesses. However, according to the compressional rheometry measurements [9], the relaxation-time of this gelatin mixture at 20Hz, is 0.6ms while that of the inclusion is approximately 3.6ms. The co- efficients of normal viscosity for 1D projection of the wave can be calculated by multiplying the relaxation-time and the elasticity. According to (5.15), this value should be divided by 3 to yield the viscosity of the material which can be calculated to be 6Pas and 36Pas for the background and inclusion, respectively. The phase of the experimental data is shown along with the re- sponses of three different simulated models constructed based on the known geometry. As expected in this case, since the thickness of the phantom is 143 Chapter 5. Longitudinal Wave Propagation smaller than the other two dimensions, a plane-stress model may lead to smaller error compared to the other two models. 5.6 Discussion The speed of propagation of the longitudinal wave is highly dependent on the boundary conditions as well as the material properties. In incompress- ible materials, when the boundary conditions force the wave to propagate only through a change of volume, a high speed (cp) p-wave occurs. The p-wave, traveling as fast as ultrasound in the medium, cannot be tracked by means of ultrasound, hence compressional transient elastography has re- ceived very limited attention in the literature. However, it is shown here that with proper boundary conditions and excitation, it is feasible to induce a longitudinally propagating wave in the medium that has a much lower speed than the p-wave. Specifically, it is shown that if a block of material can freely preserve its volume, having a large uniform exciter on the sur- face can approximate the semi-infinite thin rod situation in which an elastic wave travels at only a few meters per second (c0). On the other hand, if the large exciter is replaced by a point-source indenter, the speed can be nearly reduced to the shear wave speed in the medium (cps). One can expect that using a mid-size exciter can produce other velocities between cps and c0 while tuning the boundary conditions and the geometry can yield velocities between c0 and cp. In such complex cases, the velocity may not necessarily be constant throughout the medium. In a block of material, when the sides are confined, the volume change would be inevitable and therefore, it is only the p-wave that conveys the information. Due to near incompressibility of the material, the p-wave generated in the medium will be much faster than the elastic wave in an unbounded medium. Finite element analysis was performed to validate the theoretical findings and explore the the effect of the boundary conditions. To achieve stable finite element solutions, the Poisson’s ratio is chosen smaller than the actual value. Although a closer value to 0.5 is more realistic for incompressible soft tissues, it causes locking of the elements and instability in solving equation (5.44). Locking or hourgalssing of the elements make them behave in a stiffer fashion. The problem of locking can be reduced by numerical integration instead of closed-form integration [20]; however, this may introduce inaccuracy in the final solution. A remedy to this trade-off is to assume near incompressibility, thus ν = 0.495. As a result, locking would not happen, however the speed of the p-wave would be smaller in the simula- 144 Chapter 5. Longitudinal Wave Propagation 0 10 20 30 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 Depth (mm) Ph as e (ra d) Experiment Plane−strain Plane−stress 1D MSD Inclusion Location Figure 5.10: Phase of the displacements at 20Hz in a phantom experiment versus responses of MSD, plane-stress and plane-strain models. The phan- tom was constructed with gelatin mixture with a viscous inclusion in the middle. 145 Chapter 5. Longitudinal Wave Propagation tions compared to the actual speed of ultrasound in incompressible materials where Poison’s ratio is closer to 0.5. A closer value of Poisson’s ratio to 0.5 results in higher values of λ and cp. For example, having E = 10kPa and ν = 0.5 − 7× 10−7 gives cp = 1540m/s. The peak of the transfer functions has been used to calculate the wave speed. At the resonance frequency of a purely elastic phantom, all the points vibrate at the same phase. The first mode corresponds to the frequency that a full wavelength is trapped in the finite medium, therefore the wavelength of the wave would be twice as large as the axial dimension of the phantom. At such a frequency, the middle of the phantom experiences the highest amplitude of vibration which corresponds to a peak in the corresponding transfer function magnitude as shown in Fig. 5.4(a). The phase also drops at the natural frequency as seen in Fig. 5.4(b). The significant phase decrease in the simulations is a result of using a linear model, while real tissue often acts in nonlinear ways. As an example, human tissue and gelatin-based phantoms are known to have a power-law viscous behavior, due to which, the viscosity and relaxation-time drop as the frequency increases [9,26–28]. This nonlinearity causes the phase to have smaller change at high frequencies, however a linear FEM model is unable to account for such effects. As shown in Fig. 5.5, having a lower viscosity gives rise to various resonances in the system at high frequencies. Therefore, generally one may see several peaks in the transfer function, given a wide-band excitation. It is known that when a uniform quasi-static compression is applied to a block of material from one side while the other side is bounded, the axial displacements attenuate almost linearly. However, in the case of a point-source exciter, the intensity falls off similar to an ultrasound beam as the square of the distance from the source [29]. Therefore, the quasi- static displacement amplitude for the case of a point-source is expected to be less than that of a large exciter. This explains the smaller transfer function magnitude for a point-source compared to the other two cases in Fig. 5.4(a) and the non-uniform spacing of the nodal displacement curves in Fig. 5.6(e) while Figs. 5.6(a) and 5.6(c) exhibit uniform variation between the displacements at the tail of the curves. It has been shown that a longitudinal wave may travel at a slow speed between c0 and cps in finite medium, depending on the boundary conditions. As seen in Fig. 5.6(b) and 5.6(f), for a finite block of material a low speed elastic wave is not the only wave present in the medium. Upon exerting an axial compression to a finite block, a p-wave is first generated in the medium which travels at the high speed of cp. At this instant, the wave has not reached the boundaries and therefore, the boundary conditions and 146 Chapter 5. Longitudinal Wave Propagation geometry are not known. Once such information is collected by the p-wave, depending on the boundary conditions, a slower wave starts to propagate that is influenced by the material properties as well as the boundary 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 propagates solely through volumetric change of the medium at a speed close to cp. As noted before, in real situations where Poisson’s ratio is closer to 0.5, the speed of the p-wave is cp ≈ 1540m/s, while the speed of the slower wave is not sensitive to the Poisson’s ratio. However, in a bounded finite medium which may also represent the case of a point-source in an infinite medium, the p-wave is the only type of longitudinal wave that can be seen in Fig. 5.6(c). Note that in this case, the resonance of the p-wave lasts for a longer time compared to the p-wave resonance in Figs. 5.6(a) and 5.6(e), since the boundary conditions assert the presence of this type of wave in the medium. Although the 3D finite element method is considered as the gold stan- dard in analyzing the mechanics of the continuum, often due to extensive computational requirements a simplified 2D or 1D model is utilized to study the tissue deformations. A 2D FEM model can be formulated by making either the plane-strain or the plane-stress assumption. A plane-strain as- sumption is ideal for the case when the region has infinite thickness while a plane-stress case is ideal when the thickness is zero. Also, the viscoelastic wave equation can be approximated in 1D and discretized by a series of mass-spring-damper elements [19]. In order for the viscoelastic behavior of the 1D model to replicate that of the 3D model closely, the normal viscosity coefficient in 1D should be three times the value of the shear viscosity as expressed by equation (5.15). To what extent can the wave equation and deformations obtained from simplified models be valid, is a question that was explored in Figs. 5.8 and 5.9. For the sample geometry depicted in Fig. 5.3 with homogeneous mechan- ical properties, the transfer functions in Fig. 5.8 show that the spectral and temporal behaviors of all of the four models are comparable for frequencies as high as the natural frequency of the system. One can say that in this exam- ple, the plane-stress or even 1D model is more accurate than the plane-strain case to account for the high frequency response of the nodal displacements. Note that the transfer functions in Fig. 5.8(a) are only for the center point of the system, the displacement of which seems to have small dependency on the type of the model in Fig. 5.9(b) and 5.9(c). A low-frequency analy- sis in Fig.5.9 shows that for the same geometry and boundary conditions, 147 Chapter 5. Longitudinal Wave Propagation a plane-strain model can describe the nodal displacements and strains in a more accurate way overall. At the inclusion boundaries, the plane-stress and 1D models exhibit sharp edges while smooth transitions produced by the 3D model can be also produced by having a plane-strain assumption. However, in a phantom experiment where the thickness was smaller than the other two dimensions, plane-stress outperformed the other two models in Fig. 5.10. 5.7 Conclusion It is feasible to use longitudinal waves in order to identify the mechanical properties in soft tissues. While the p-wave cannot be tracked with ul- trasound, low-speed longitudinal waves can be induced by controlling the excitation and boundary conditions. In this case, first a fast primary wave is generated in the medium and after it reaches the boundaries, a slow sec- ondary wave will emerge. The wave equation can be discretized by a dynamic finite element model or it can be simplified to 1D by a series of mass-spring damper elements. However, the 1D wave equation fails to account for some critical phenomena in the medium. The viscosity coefficient that should be incorporated in the 1D equation is also different from the original viscosity that is considered in 2D and 3D models, being almost three times larger in the case of incompressible materials with zero bulk viscosity. 148 References [1] S. Catheline, J. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelka- ram, and J. Culioli, “Measurement of viscoelastic properties of homo- geneous soft solid using transient elastography: An inverse problem ap- proach,” Journal of the Acoustical Society of America, vol. 116, no. 6, pp. 3734–41, 2004. [2] M. Zhang, B. Castaneda, Z. Wu, P. Nigwekar, J. V. Joseph, D. J. Rubens, and K. J. Parker, “Congruence of imaging estimators and mechanical measurements of viscoelastic properties of soft tissues”, Ultrasound in Medicine and Biology, vol. 33, no. 10, pp. 1617–31, 2007. [3] M. Fatemi and J. F. Greenleaf, “Vibro-acoustography: An imaging modality based on ultrasound-stimulated acoustic emission,” Proceed- ings of the National Academy of Science, U.S.A., vol. 96, no. 12, pp. 6603–08, June 1999. [4] J. Bercoff, M. Tanter, and M. Fink, “Supersonic shear imaging: A new technique for soft tissue elasticity mapping,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 51, no. 4, pp. 396–409, 2004. [5] A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov, “Shear wave elasticity imaging: a new ultrasonic tech- nology of medical diagnostics,” Ultrasound in Medicine and Biology, vol. 24, no. 9, pp. 1419–35, 1998. [6] K. Nightingale, S. McAleavey and G. Trahey, “Shear-wave generation using acoustic radiation force: in vivo and ex vivo results”, Ultrasound in Medicine and Biology, vol. 29, no. 12, pp. 1715–23, 2003. [7] J. McLaughlin, and D. Renzi, “Shear wave speed recovery in transient elastography and supersonic imaging using propagating fronts”, Inverse Problems, vol. 22, no. 2, pp. 681–706, 2006. 149 Chapter 5. Longitudinal Wave Propagation [8] R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink, “Viscoelastic shear properties of in vivo breast lesions measured by MR elastography,” Magnetic Resonance Imaging, vol. 23, no. 2, pp. 159–65, 2005. [9] H. Eskandari, S. Salcudean, and R. Rohling, “Viscoelastic parameter estimation based on spectral analysis,” IEEE Transactions on Ultrason- ics, Ferroelectrics and Frequency Control, vol. 55, no. 7, pp. 1611–25, July 2008. [10] A. Baghani, H. Eskandari, S. Salcudean, and R. Rohling, “Measure- ment of Viscoelastic Properties of Tissue Mimicking Material Using Longitudinal Wave Excitation,” Submitted to IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Nov. 2008. [11] M. F. Insana, C. Pellot-Barakat, M. Sridhar, and K. K. Lindfors, “Vis- coelastic imaging of breast tumor microenvironment with ultrasound,” Journal of Mammary Gland Biology and Neoplasia, vol. 9, no. 4, pp. 393–404, Oct. 2004. [12] M. Insana, J. Liu, M. Sridhar and C. Pellot-Barakat, “Ultrasonic me- chanical relaxation imaging and the material science of breast cancer”, Proceedings of the IEEE Ultrasonics Symposium, vol. 2, pp. 739–42, 2005. [13] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, “Time-resolved 2D pulsed elastography: experiments on tissue-equivalent phantoms and breast in vivo”, Proceedings of SPIE, vol. 4325, pp. 120–26, 2001. [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 Ultrasonics Ferroelectrics and Frequency Control vol. 49, no. 4, pp. 436–46, 2002. [15] L. Sandrin, D. Cassereau, and M. Fink, “The role of the coupling term in transient elastography”, J. Acoust Society America, vol. 115, no. 1, pp. 73–83, 2004. [16] J. H. C. Thompson, “On the theory of visco-elasticity: A thermody- namical treatment of visco-elasticity, and some problems of the vibra- tions of visco-elastic solids”, Philosophical Transactions A, vol. 231, pp. 339–407, 1933. 150 Chapter 5. Longitudinal Wave Propagation [17] L. E. Malvern, Introduction to the mechanics of a continuous medium, 1st ed. New Jersey: Prentice-Hall Inc., 1969. [18] H. Kolsky, Stress Waves in Solids. Dover Publications Inc., 1963. [19] E. Turgay, S. E. Salcudean, and R. Rohling, “Identifying the mechan- ical properties of tissue by ultrasound strain imaging,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 221–35, 2006. [20] R. D. Cook, D. S. Malkus, and M. E. Palesha, Concepts and Applica- tions of Finite Element Analysis, 3rd ed. New York: John Wiley & Sons, Inc., 1989. [21] H. Eskandari, S. Salcudean, R. Rohling and J. Ohayon, “Viscoelastic characterization of soft tissue from dynamic finite element models,” Physics in Medicine and Biology, vol. 53, no. 22, pp. 6569-90, Nov. 2008. [22] M. Bro-Nielsen, “Finite element modeling in surgery simulation,” Pro- ceedings of the IEEE, vol. 86, no. 3, pp. 490–503, March 1998. [23] D. Fu, S. F. Levinson, S. M. Gracewski, and K. J. Parker, “Non-invasive quantitative reconstruction of tissue elasticity using an iterative forward approach,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1495– 509, Aug 2000. [24] R. Zahiri-Azar and S. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, Oct. 2006. [25] L. Ljung, System identification : theory for the user, Prentice Hall PTR, NJ, 1999. [26] M. Kiss, T. Varghese, and T. Hall, “Viscoelastic characterization of in vitro canine tissue,” Physics in Medicine and Biology, vol. 49, no. 18, pp. 4207–18, Sept. 2004. [27] C. Joly-Duhamel, D. Hellio, and M. Djabourov, “All gelatin networks: 1. Biodiversity and physical chemistry,” Langmuir, vol. 18, no. 19, pp. 7208–17, 2002. [28] N. Phan-Thien, S. Nasseri, and L. Bilston, “Oscillatory squeezing flow of a biological material,” Rheologica Acta, vol. 39, no. 4, pp. 409–17, 2000. 151 Chapter 5. Longitudinal Wave Propagation [29] R. S. C. Cobbold,, Foundations of Biomedical Ultrasound, Oxford Uni- versity Press, NY, 2007. 152 Chapter 6 Conclusions and Future Research In Chapters 2, 3 and 4 novel algorithms were proposed to implement a complete vibro-elastography system. All the proposed methodologies were tested on simulated data as well as tissue-mimicking phantoms. In Chapter 5 the theoretical framework for vibro-elastography and any other compres- sional viscoelasticity imaging set-up was presented and feasibility of such techniques was proved. In this chapter, a full analysis of the materials in the previous chapters is presented in light of current research in the field. The strengths and weaknesses of the approaches undertaken in this thesis are described and the the contributions of the thesis are delineated. The future work and pathways are discussed in the end. At this point, it can be deduced that the hypothesis of this thesis as stated in Chapter 1 is confirmed through numerical simulations and exper- iments on tissue-mimicking phantoms. Hence, it is feasible to estimate the viscoelastic properties of soft materials from the dynamic tissue response to a transient or steady-state longitudinal or shear excitation. 6.1 Contributions 6.1.1 Strain Imaging With the development of a strain imaging algorithm proposed in Chapter 2, objective 1 of the thesis was met and it was made possible to explore ob- jective 2 and objective 3 of the thesis through phantom experiments. Any ultrasound-based elastography system entails a motion estimation subsys- tem. As the first objective of the thesis, a new strain imaging algorithm was developed in Chapter 2 to find and match the peaks in pre- and post- compression RF A-lines. A significant novelty of this algorithm is the use of the wavelet transform to identify the true peaks of the signals. With the specific choice of the mother wavelet function, applying the wavelet trans- form on an RF signal results in low-pass filtering and then differentiating 153 Chapter 6. Conclusions and Future Research the signal. The local compression of the RF signals was associated with the displacement of the peaks. The peaks of the pre- and post-compression RF signals were then identified through locating the zero-crossings of the wavelet transforms. In comparison to the zero-crossing algorithm by Srinivasan and Ophir [1], the wavelet filtering of PSA results in better noise rejection and lower pos- sibility of detecting the false peaks. Compared to the conventional cross- correlation method, the proposed algorithm has a better signal-to-noise ra- tio, level of error, and speed of processing. Basically, the PSA method uses a feature of narrow-band RF signals which is the high number of peaks per each line, while cross-correlation techniques utilize a general definition of similarity which can be common to most types of signals. The PSA algo- rithm, being particularly designed and tuned for ultrasound RF signal mo- tion estimation, is therefore expected to result in high sensitivity and SNR. The relatively high efficiency of PSA is due to its simplicity, which only in- volves three steps: convolution, zero-crossing detection and peak matching. The wavelet transform can be implemented in the form of a convolution with the mother wavelet function. Detection of the zero-crossings is a fast process, however subsample allocation of the peaks or wavelet transform zero-crossings involves an extra step which can also be a fast algorithm. A sophisticated technique was proposed for the peak matching algorithm. The peak matching step is critical and should be implemented accurately. If one of the post-compression peaks is incorrectly associated to a peak in the pre-compression signal, there is a high chance that the rest of the peaks will be incorrectly matched and the algorithm fails. From this aspect, the cross-correlation technique can be implemented in a more efficient way by considering prior estimates from previous blocks [2]. This way, the cross- correlation technique can achieve very high frame-rates and possibly perform faster, but still less accurate than the PSA. PSA has the potential for increasing the resolution to the level of inter- peak distance by reducing the window length while maintaining the overlap value. Therefore, the resolution of PSA may be increased to the degree of the wavelength of the ultrasound beam. If strain is sufficiently small and the peaks are found and matched correctly, the resolution can be reduced to the level of the ultrasound wavelength for the proposed method, while such a resolution is not achievable with most other methods. This increase in resolution will only be achieved at the cost of lowering the signal-to-noise ratio. 154 Chapter 6. Conclusions and Future Research 6.1.2 Spectral Viscoelastic Parameter Estimation The spectral algorithm to estimate the viscoelastic properties of soft materi- als proposed in Chapter 3 satisfied objective 2 of the thesis, while objective 5 was also met through the phantom construction proposed in that chapter. A number of original contributions were made in this work. The first contri- bution was using the phase of the transfer functions or spectra to calculate the relaxation-time. The idea originated from the concept of the transfer functions as suggested by Turgay et. al [3]. The algorithm builds upon the discretization that was used in [3]. How- ever, in Chapter 3 it was proposed that the model can be modified for low-frequency measurements. The relaxation-time can be estimated from the phase of the strains or displacements. The phase from either a single- frequency mechanical excitation or the phase from a multiple-frequency transfer function approach can be used to estimate the relaxation-time. Once the relaxation-time is estimated, one possible step is to fit the dataset of relaxation-time and frequency to a power-law relationship, deriving power- law parameters that can be estimated for further tissue characterization. Power-law represents a nonlinear viscous effect in the frequency domain, which can be modeled as a Kelvin-Voigt model with fractional derivatives in the time domain [4] and has been reportedly observed in soft biological materials [5]. The proposed method provides a novel way of calculating the relaxation-time from the phase of the strain or displacement spectra, with the possibility of estimating the power-law behavior of the relaxation-time and/or viscosity as an additional feature for tissue characterization. The relaxation-time is known to be of particular significance when dif- ferent types of material need to be distinguished. For example, a soft gelatin phantom with low concentration of gelatin powder was seen to have the same relaxation-time as a hard gel of the same kind with relatively high gelatin powder concentration. In order to test the derived methods, a PVA-based sponge-like material which was observed to have a high damping was embed- ded in a gelatin phantom to model an inclusion with high relaxation-time and viscosity. Using such a technique to fabricate tissue-mimicking materi- als for elastography experiments was another original contribution that was made in that work. An experimental apparatus and software were developed to test the al- gorithms. The setup can be used both as a rheometer to characterize the overall mechanical properties of a material or as a vibro-elastography imag- ing device using an ultrasound system. In addition to the numerical sim- ulations, the accuracy of the estimation methods were validated with this 155 Chapter 6. Conclusions and Future Research compressional rheometer. The rhometry tests confirmed the elasticity and relaxation-time contrasts as estimated with the proposed spectral method of ultrasound vibro-elastography. A drawback of the method proposed in this work is the use of a simple 1D model to account for the tissue deformation. As a result, the effects of the boundary conditions and excitation profile are neglected. Application of this algorithm to in vivo situations requires that a flat exciter applies a symmetric compression. The boundary conditions, including the neighbor- ing inclusions should also be symmetric to minimize lateral and out-of-plane motion. Also, bulging should be permissible within the compressed region, so that a low-speed elastic was can propagate in the tissue and tracked by ultrasound and the a significant phase can be detected between local displacements and also between axial strains. However, the proposed esti- mation algorithm has the potential to be efficiently implemented to image the viscoelastic properties of soft tissue in real-time. 6.1.3 Finite Element Simulations In order to handle the effects of boundary conditions on three dimensional deformations of the tissue, a dynamic finite element model with a novel de- finition of the damping matrix was proposed in Chapter 4. Finite element methods have been widely explored in the context of static elastography problems [6,7]. With the increasing importance of viscosity and damping in viscoelasticity imaging, dynamic finite element models are gaining more in- terest in recent years [8–10]. The damping models that are presently used for modeling the soft tissue are generally adopted from the structural mechanics where a linear combination of mass and stiffness matrices is used to account for the damping in the structures [10, 11]. In Chapter 4, a new method to form the damping matrix for finite element analysis has been proposed where the effective damping is assumed to originate from the viscosity of the medium. Thus, the viscous deformations would follows a Voigt’s viscoelastic deformation model of soft tissue. The model assumes known lumped masses at the nodes, and comprises two vectors of elasticity and viscosity para- meters that depend on the material elasticity and viscosity distributions, respectively. 6.1.4 Finite Element Inverse Problems Solving the inverse problem of viscoelasticity in Chapter 4 satisfied objec- tive 3 of the thesis. An algorithm was devised to estimate the elasticity 156 Chapter 6. Conclusions and Future Research and viscosity in the medium using the complex displacement field in a finite element platform. An inverse algorithm was formulated that searches for the optimal values of the Youngs modulus and viscosity that yield the least deviation from the observed displacements inside the medium. Using this deformation model and the observed dynamic data of a harmonic excitation test, the inverse problem was solved to reconstruct the viscosity and elas- ticity in the medium by using a Gauss-Newton based approach. As in other inverse problems, previous knowledge of the parameters on the boundaries of the medium are necessary to assure uniqueness and convergence and to obtain an accurate map of the viscoelastic properties. If 3D data is not available, the problem can be simplified to 2D by as- suming plane-strain, plane-stress or any other type of 2D deformation state. If the lateral displacements are not available or have a a poor quality, the method can be formulated in terms of the axial data only. However, proper and reasonable assumptions for the lateral displacements or forces on the boundaries of the ROI are necessary to converge to the correct distribution of the parameters. To achieve this, it is preferable to have lateral displace- ments estimated at least on the boundaries of the ROI. In in vivo exams, a priori knowledge of the parameters at some point inside the ROI is required to converge to the correct values of elasticity and viscosity. An inaccurate initial condition assumption may result in the scaling of the elasticity, while the viscosity estimate may become biased and inaccurate. 6.1.5 Justifying the Feasibility of Compressional Elastography The results of the work in Chapter 5 helped meet objective 4 of the thesis. While applying shear wave was assumed to be the only feasible method of dynamic tissue characterization with ultrasound or MRI, Chapters 3 and 4 indicated that longitudinal excitation can also be used to estimate the mechanical properties of tissue. The general assumption in the literature is that longitudinal wave propagates at the speed of ultrasound, in which case it would be impossible to image the deformations with ultrasound mo- tion estimation techniques. The theoretical foundation for the feasibility of the longitudinal wave elastography, upon which the concept of vibro- elastography is established is presented in Chapter 5. It is proved that the speed of propagation is highly dependent on the boundary conditions. Therefore, the longitudinal wave speed can have any value between that of s- wave (shear wave) and that of p-wave (primary wave). In a recent work [12], we have tracked an externally induced longitudinal wave in incompressible 157 Chapter 6. Conclusions and Future Research tissue-mimicking materials with ultrasound in order to estimate the Young’s modulus and viscosity from the standing wave patterns at low-frequencies (< 150Hz). The method is based on fitting the appropriate viscoelastic pa- rameters to the longitudinal wave model in which the wave speed is predicted to be √ E/ρ. 6.1.6 Analysis of the Longitudinal Wave Propagation Finite element simulations in Chapter 5 confirmed the theoretical deriva- tions of the wave speed for different boundary conditions. A 3D FEM model was used for analysis and validation purposes, while the error produced by simplified 2D FEM or even 1D models were explored. It was seen that a simplified model can be reasonable in some aspects in modeling the defor- mations, while such simplifications may result in drastic errors or biases in other respects, especially for non-homogeneous environments. The prop- agation of low speed elastic wave and high speed p-wave were studied in simulations, resulting in a better understanding of the propagation of the longitudinal waves in a homogeneous medium. 6.2 Future Work The algorithms presented in this thesis yield novel methods for viscoelastic parameter estimation in soft human tissue. If implemented in real-time, they can present a useful tool for in vivo clinical monitoring of cancer or other malignancies. As imaging the mechanical properties of soft tissue becomes more mainstream and at the same time more data are reported in the literature of the viscoelastic parameters of human organs, the potential for developing a new modality of medical imaging in the form of a new gold standard for cancer detection or real-time, low-cost tissue characterization increases significantly. The methodology and ideas that were developed in this thesis can directly enter this competing market. A number of interesting and relevant areas of future work which can be pursued from here are as follows: • Two dimensional strain imaging is possible by generalizing the idea of peak search algorithm to the 2D RF images. The peaks can be located with a subsample resolution in 2D and a sophisticated peak matching technique can be devised for more accurate estimation of the axial and lateral displacements. 158 Chapter 6. Conclusions and Future Research • If reliable transfer functions between displacements or strains can be obtained at higher frequencies, more sophisticated parametric models can be analyzed for soft tissue. Some parameters may have direct physical interpretations such as elasticity, viscosity and density, while other parameters may correspond to more subtle behaviors of human soft tissues. • Accurate two dimensional modeling of the deformations can replace standard 2D assumptions of plane-strain or plane-stress. Depend- ing on the geometry of the problem, a 2D deformation model can be formulated that describes the deformations more accurately than the standard 2D models. This way, in the FEM-based inverse problem of viscoelasticity the artifacts that result from inaccurate modeling can be reduced. • With the proposed FEM formulation, other solution strategies to the inverse problem can be investigated including non-iterative approaches. Also, the effect of lateral motion, boundary conditions and out-of-plane inhomogeneity in in vivo situations can be explored on the result of the inverse problem. • Implementation of the proposed reconstruction algorithms in real-time enables one to clinically investigate the importance of viscosity and elasticity in soft tissue characterization. For this reason special exciters can be built to exert low-amplitude, single-frequency or wide-band displacements to the surface of the body to image different organs such as breast, liver, kidney or prostate. • The developed techniques can be incorporated with an ARFI system that exerts a localized force to the tissue, internally. The wave pat- tern induced by ARFI can be investigated to see whether a slow lon- gitudinal wave can be observed. If feasible, the estimation methods developed in this thesis can be directly applied to the resulting dis- placements of an ARFI excitation. 159 References [1] S. Srinivasan and J. Ophir, “A zero-crossing strain estimator for elas- tography,” Ultrasound in Medicine and Biology, vol. 29, no. 2, pp. 227– 38, Feb 2003. [2] R. Zahiri-Azar and S. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, Oct. 2006. [3] E. Turgay, S. E. Salcudean, and R. Rohling, “Identifying the mechan- ical properties of tissue by ultrasound strain imaging,” Ultrasound in Medicine and Biology, vol. 32, no. 2, pp. 221–35, 2006. [4] M. Kiss, T. Varghese, and T. Hall, “Viscoelastic characterization of in vitro canine tissue,” Physics in Medicine and Biology, vol. 49, no. 18, pp. 4207–18, Sept. 2004. [5] N. Phan-Thien, S. Nasseri, and L. Bilston, “Oscillatory squeezing flow of a biological material,” Rheologica Acta, vol. 39, no. 4, pp. 409–17, 2000. [6] F. Kallel and M. Bertrand, “Tissue elasticity reconstruction using linear perturbation method,” IEEE Transactions on Medical Imaging, vol. 15, no. 3, pp. 299–313, Jun 1996. [7] Y. Zhu, T. J. Hall, and J. Jiang, “A finite-element approach for young’s modulus reconstruction,” IEEE Transactions on Medical Imag- ing, vol. 22, no. 7, pp. 890–901, Sep. 2003. [8] D. Fu, S. F. Levinson, S. M. Gracewski, and K. J. Parker, “Non-invasive quantitative reconstruction of tissue elasticity using an iterative forward approach,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1495– 509, Aug 2000. 160 Chapter 6. Conclusions and Future Research [9] K. Nightingale, M. Palmeri and G. Trahey, “Analysis of contrast in images generated with transient acoustic radiation force,” Ultrasound in Medicine and Biology, vol. 32, no. 1, pp. 61-72, 2006. [10] M. D. J. McGarry and E. E. W. Van Houten, “Use of a Rayleigh damp- ing model in elastography,” Medical and Biological Engineering and Computing, vol. 46, no. 8, pp. 759–66, June 2008. [11] D. F. Pilkey, G. Park, and D. J. Inman, “Damping matrix identification and experimental verification,” in Proceedings of SPIE, vol. 3672, 1999, pp. 350–357. [12] A. Baghani, H. Eskandari, S. E. Salcudean and R. Rohling, “Measure- ment of Viscoelastic Properties of Tissue Mimicking Material Using Longitudinal Wave Excitation”, Submitted to IEEE Trans. on Ultra- sonics, Ferroelectrics and Frequency Control, Nov. 2008. 161 Appendix A Wavelet Transforms and Peak Detection A complex-valued function ψ(x) of a real argument x is said to be a wavelet function or mother wavelet if its Fourier transform Ψ(ω) satisfies, ∫ +∞ 0 |Ψ(ω)|2 ω dω = ∫ 0 −∞ |ψ(ω)|2 ω dω < +∞ . (A.1) The wavelet transform (WT) of the function f(x) ∈ L2(R) at scale s is defined by the following convolution [23]: Wsf(x) = f(x) ∗ ψs(x) , (A.2) where the function ψs(x) is the dilation of the wavelet function ψ(x) and is defined as, ψs(x) = 1 s ψ(x/s) , (A.3) From (A.2) it can be seen that the continuous wavelet transform (CWT) of a function at a certain scale involves filtering of the original function. From (A.3), as s increases, the dilation of the wavelet function becomes wider in the time domain. Thus, the CWT of a signal shows the detailed or higher frequency components of that signal in its lower scales and low frequencies can be observed in the higher scales. If the scales are discretized such that s = 2j(j ∈ Z), the resulting WT is called the dyadic wavelet transform (DWT). As Mallat et. al showed in [22], if the wavelet function is the first or higher order derivative of a smoothing function, the singularities and zero- crossings of the wavelet transform would correspond to the peaks, zeros or inflection points of the original signal. A smoothing function θ(x) is defined as a function that converges to zero at infinity and its integral is equal to one. Such a function can smooth the fluctuations of the signal at different scales. To show how the CWT of a function f(x) behaves in this case, one 162 Appendix A. Wavelet Transforms and Peak Detection 0 0.5 1 1.5 2 x 107Frequency (Hz) N or m al iz ed A m pl itu de j=0 j=1j=2j=3 (a) 0 0.5 1 1.5 2 x 107Frequency (Hz) Am pl itu de (b) Figure A.1: (a) Frequency response of the equivalent filtering effect of the mother wavelet shown in Fig. 6(c) when dilated on the dyadic scales. The amplitude is normalized to get the best clarification. (b) Fourier Trans- form of an ultrasound RF line acquired using a 2.5MHz probe with 40MHz sampling rate. 163 Appendix A. Wavelet Transforms and Peak Detection can define a wavelet function ψa(x) as the first derivative of a smoothing function θ(x): ψa(x) = dθ(x) dx . (A.4) It can be shown that the zero integral condition for the wavelet functions is satisfied by this choice of the mother wavelet. The CWT of f(x) using this wavelet function is: W as f(x) = f(x) ∗ ψas (x) = s d dx (f ∗ θs)(x) , (A.5) where ∗ is the convolution operator and the index a indicates that the wavelet transform has been calculated using ψa(x) as the mother wavelet. Equation (A.5) shows that the local extrema of the smoothed signal f ∗θ(x) correspond to the zero-crossings in W as f(x), whereas the inflection points in the smoothed signal will result in the extrema in W as f(x) , as has been illustrated in Fig. 2.2. The maxima and the minima of a continuous wavelet transform are generally referred to as modulus maxima. From equation (A.5), the effect of the CWT with the choice of wavelet function shown in Fig. 2.2(c), on a signal is to filter it in a certain frequency range as defined by the dilation of the wavelet function, and then taking the derivative. The frequency response of the smoothing function when dilated on different dyadic scales is shown in Fig. A.1(a). Also, a typical spectrum of an ultrasound RF A-line is shown in Fig. A.1(a). The signal shown in this example results from a 2.5MHz transducer, scanning through a homogeneous phantom. This figure shows that for retaining most of the energy of the signal and still filtering the high frequency noise and low frequency interferences, the best choice would be scale 22 that has the most overlap with the spectrum of the RF signals. By knowing the spectrum of the ultrasound RF signal, a suitable scale can be selected in order to get the best filtering effect. For high frequency transducers, even non-dyadic scales that match the spectrum of the RF signal can be selected as well. 164 Appendix B Regularity of the Wavelet Transform A function f(x) is said to be uniformly Lipschitz of order α, 0 < α < 1, over an open interval ]a, b[, if and only if there exists a constant K > 0 such that for any (x0, x1) ∈]a, b[2, |f(x0)− f(x1)| ≤ K |x0 − x1|α . (B.1) The Lipschitz uniform regularity of f(x) is the upper bound α0 of all α that satisfy the above inequality [27]. With a similar definition, it can be shown that the above equation can be expressed in terms of the wavelet transform of the signal as follows: |W2jf(x)| ≤ K(2j)α . (B.2) Higher regularity implies that the signal has more singular characteristics at that point. Li et al. [9] defined the regularity of a signal at a certain point as the ratio of the WT values of that signal at two consecutive scales. Likewise, we define the regularity of a point x0 in the function f(x) as: α = log2(W2j+1f(x0))− log2(W2jf(x0)) . (B.3) Due to high frequency interference, there might be several false peaks in the RF signal. In order to find the true peaks, one can use the concept of regularity. Based on the energy spectrum of the RF signal, a true peak can be defined as one with most of its energy in scale 21 while a false peak is one with more energy in scale 20. This means that the tangent of the wavelet transform of a true peak is higher in scale 20 than 21 while for a false peak, the opposite applies. Therefore, letting j = 0, if α in (B.3) becomes greater than zero for both modulus maxima of a certain peak, that would be a true peak location; otherwise it would be a false peak and should be eliminated from further processing. In dynamic elastography and cases in which speed is most important, this step can be ignored. 165 Appendix C Transfer Function Parameterization The constitutive equation of a Voigt element can be obtained in the fre- quency domain by taking the Fourier transform of eqn (3.1): (E + jωη)ǫ̂(ω) = σ̂(ω) , (C.1) where ǫ(ω) is the Fourier transform of the strain at a specific location in tissue (any node in Fig. 3.1) and σ(ω) is the Fourier transform of the stress signal. E and η are the local values of the elasticity and viscosity. The transfer function between the strain and the stress is defined as: Hǫσ(ω) = Pǫσ(ω) Pσσ(ω) , (C.2) where Pǫσ(ω) is the cross-power spectral density between ǫ(t) and σ(t) and Pσσ(ω) is the power spectral density for σ(t). If stress has been recorded, eqn (C.1) can be expressed as: (E + jωη) = 1 Hǫσ(ω) . (C.3) Otherwise, the stress can be removed from eqn (C.1) as follows: (E + jωη)ǫ(ω) = (E0 + jωη0)ǫ0(ω) , (C.4) where E0 and η0 are the parameters at any arbitrary location (along the 1D line) and ǫ0 is the strain at that point. By taking the transfer function between the strains at the two points, the following results: (E + jωη) (E0 + jωη0) = 1 Hǫǫ0(ω) . (C.5) Different characteristics of the transfer functions, such as the magnitude 166 Appendix C. Transfer Function Parameterization and the phase depend on the tissue properties. Identification of these me- chanical properties from the transfer functions results in good estimates, since firstly, the amount of temporal averaging to calculate the power spec- tral density functions can be increased arbitrarily and secondly, the transfer function can be explored at several frequencies to extract a specific feature. Instead of using the transfer function between the stress and the strain, one can use the applied force and the nodal displacements to calculate the transfer functions. From eqn (3.2), the following results: (ki + jωbi) = 1 H if (ω)−H i+1f (ω) , (C.6) where H if (ω) denotes the transfer function between the force and the dis- placement at node i. The stiffness and the relaxation-time can be calculated here, similarly to the case in which strains and stress are used to compute the transfer functions. 167 Appendix D Elasticity Estimation from the Transfer Functions It has been shown in [11] that using eqn (C.3), the following estimate for the elasticity can be obtained: Ê = 1∣∣H̄ǫσ(ω)∣∣ , (D.1) where Ê is the estimated elastic modulus and H̄ǫσ(ω) is the low-frequency asymptote of the transfer function. The term jωη in eqn (C.3) is ignored due to small frequency approximation and the fact that ωη << E. If the stress is not measured, eqn (C.5) can be used to get the following relative estimate for the elasticity: Ê = E0∣∣H̄ǫǫ0(ω)∣∣ . (D.2) The estimates should be normalized with respect to E0 which is unknown. Therefore, only the profile of the relative elasticity can be obtained. 168 Appendix E Simulating a Mass-Spring-Damper Network The dynamic motion of the tissue can be described by a series of mass- spring-damper elements as depicted in Fig. 3.1. In such a network, the tissue motion is governed by the following equation [11]: KX(t) +BẊ(t) +MẌ(t) = F (t) , (E.1) where X(t) is the matrix that contains the displacements of all the nodes versus time and F (t) is the external force. The matrices K, B and M are the stiffness, damping and mass matrices, defined in [11]. Therefore, a state-space model can be generated to compute the nodal displacement vectors: [ Ẋ(t) Ẍ(t) ] = [ 0 I −M−1K −M−1B ] [ X(t) Ẋ(t) ] + [ 0 M−1 ] F (t) , X(t) = [ I 0 ] [ X(t) Ẋ(t) ] . (E.2) The output, X(t), can be numerically evaluated by direct discretization of the above model, using a fixed step-size. Therefore, the Matlab function (Mathworks Inc., Natick, MA, USA), lsim, has been used in this paper and the sampling frequency was fixed at 1KHz. 169
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- On the identification of mechanical properties of viscoelastic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
On the identification of mechanical properties of viscoelastic materials Eskandari, Hani 2009
pdf
Page Metadata
Item Metadata
Title | On the identification of mechanical properties of viscoelastic materials |
Creator |
Eskandari, Hani |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | Commonly used medical imaging techniques can render many properties of the anatomy or function, but are still limited in their ability to remotely measure tissue mechanical properties such as elasticity and viscosity. A remote and objective palpation function would help physicians in locating possible tumors or malignancies. The branch of medical imaging that characterizes tissues mechanical properties in a non-invasive manner has enjoyed increasing interest in the past two decades. The basic principle is to apply an excitation, such as tissue compression, to a region of interest and measure the resulting tissue deformation. Tissue mechanical properties can then be inferred from the observed deformation at multiple locations in the region, and the properties can be displayed as an image. If the excitation is dynamic, the deformation is considered as a motion field that varies in time and location over the region of interest. Ultrasound is particularly well suited for measuring motion fields due to its ability to image in real-time, low cost, low risk and ease-of-accessibility. The focus of this thesis is the estimation of the viscoelastic parameters such as Young's modulus, viscosity and relaxation-time. For this purpose, a motion estimation method is proposed to measure axial tissue displacements from the peak of the ultrasound radio frequency signals. The displacements can be further processed to identify the mechanical properties. Two methods were developed: the first one is based on a one dimensional Voigt's model of soft tissue and the second one is based on a finite element model. In the first method, a single frequency or wide-band excitation is applied to the tissue and the local relaxation-time is recovered from the phase difference between the strains or displacements. In this method, the elasticity can also be reconstructed from the magnitudes of the spectra. In the second approach, a novel dynamic finite element model is proposed for the incompressible soft materials. An inverse problem of viscoelasticity is solved iteratively to reconstruct the viscosity and elasticity based on a two or three dimensional model. The theoretical aspect of compressional elastography and longitudinal wave propagation is investigated. It is shown to be feasible to apply dynamic or transient compressional excitation to recover the dynamic properties of soft tissue. |
Extent | 1447107 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067075 |
URI | http://hdl.handle.net/2429/6051 |
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 |
GraduationDate | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_spring_eskandari_hani.pdf [ 1.38MB ]
- Metadata
- JSON: 24-1.0067075.json
- JSON-LD: 24-1.0067075-ld.json
- RDF/XML (Pretty): 24-1.0067075-rdf.xml
- RDF/JSON: 24-1.0067075-rdf.json
- Turtle: 24-1.0067075-turtle.txt
- N-Triples: 24-1.0067075-rdf-ntriples.txt
- Original Record: 24-1.0067075-source.json
- Full Text
- 24-1.0067075-fulltext.txt
- Citation
- 24-1.0067075.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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067075/manifest