SHEAR ELASTICITY AND SHEAR VISCOSITY IMAGING IN SOFT TISSUE
By
Yiqun Yang

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Electrical Engineering – Doctor of Philosophy
2018

ABSTRACT
SHEAR ELASTICITY AND SHEAR VISCOSITY IMAGING IN SOFT TISSUE
By
Yiqun Yang
In this thesis, a new approach is introduced that provides estimates of shear elasticity and
shear viscosity using time-domain measurements of shear waves in viscoelastic media.
Simulations of shear wave particle displacements induced by an acoustic radiation force
are accelerated significantly by a GPU. The acoustic radiation force is first calculated using
the fast near field method (FNM) and the angular spectrum approach (ASA). The shear
waves induced by the acoustic radiation force are then simulated in elastic and viscoelastic
media using Green’s functions. A parallel algorithm is developed to perform these calculations on a GPU, where the shear wave particle displacements at different observation points
are calculated in parallel. The resulting speed increase enables rapid evaluation of shear
waves at discrete points, in 2D planes, and for push beams with different spatial samplings
and for different values of the f-number (f/#). The results of these simulations show that
push beams with smaller f/# require a higher spatial sampling rate. The significant amount
of acceleration achieved by this approach suggests that shear wave simulations with the
Green’s function approach are ideally suited for high-performance GPUs.
Shear wave elasticity imaging determines the mechanical parameters of soft tissue by
analyzing measured shear waves induced by an acoustic radiation force. To estimate the
shear elasticity value, the widely used time-of-flight (TOF) method calculates the correlation between shear wave particle velocities at adjacent lateral observation points. Although
this method provides accurate estimates of the shear elasticity in purely elastic media, our
experience suggests that the TOF method consistently overestimates the shear elasticity
values in viscoelastic media because the combined effects of diffraction, attenuation, and
dispersion are not considered. To address this problem, we have developed an approach

that directly accounts for all of these effects when estimating the shear elasticity. This new
approach simulates shear wave particle velocities using a Green’s function-based approach
for the Voigt model, where the shear elasticity and viscosity values are estimated using an
optimization-based approach that compares measured shear wave particle velocities with
simulated shear wave particle velocities in the time-domain. The results are evaluated on
a point-by-point basis to generate images. There is good agreement between the simulated
and measured shear wave particle velocities, where the new approach yields much better
images of the shear elasticity and shear viscosity than the TOF method.
The new estimation approach is accelerated with an approximate viscoelastic Green’s
function model that is evaluated with shear wave data obtained from in vivo human livers.
Instead of calculating shear waves with combinations of different shear elasticities and shear
viscosities, shear waves are calculated with different shear elasticities on the GPU and then
convolved with a viscous loss model, which accelerates the calculation dramatically. The
shear elasticity and shear viscosity values are then estimated using an optimization-based
approach by minimizing the difference between measured and simulated shear wave particle
velocities. Shear elasticity and shear viscosity images are generated at every spatial point
in a two-dimensional (2D) field-of-view (FOV). The new approach is applied to measured
shear wave data obtained from in vivo human livers, and the results show that this new
approach successfully generates shear elasticity and shear viscosity images from this data.
The results also indicate that the shear elasticity values estimated with this approach are
significantly smaller than the values estimated with the conventional TOF method and that
the new approach demonstrates more consistent values for these estimates compared with
the TOF method. This experience suggests that the new method is an effective approach
for estimating the shear elasticity and the shear viscosity in liver and in other soft tissue.

Copyright by
YIQUN YANG
2018

Dedicated to Qianwei Jiang

v

ACKNOWLEDGEMENTS

This thesis is a summary of my Ph.D. work, but more importantly, it is also a milestone for
my five years journey at Michigan State University. During the past five years, I’ve received
a considerable amount of help from several people. Here, I would like to express my sincere
thanks to them.
First and foremost, I would like to thank my Ph.D. advisor, Prof. Robert McGough.
I would like to thank him for providing me this opportunity to work in the Biomedical
Ultrasonics and Electromagnetics Lab, and more specifically, on this very interesting project
on ultrasound elastography. With his guidance, mentoring, and knowledge, I was able to
gradually build up my expertise in medical ultrasound. I’ve also learned a lot from his
diligence, earnestness, and enthusiasm, which taught me what a researcher should be.
I would also like to thank my Ph.D. committee members: Dr. Seungik Baek, Dr. Brian
Feeny, and Dr. Lalita Udpa. I would also like to thank them for useful discussions and
suggestions on my research.
I would like to thank my collaborators at Mayo Clinic. I especially thank Dr. Matthew
Urban for sharing knowledge, explaining methods, sending data, and providing insightful
discussions for my ultrasound elastography project. I also want to thank Bo Qiang for the
great collaborations on the finite element and finite difference simulations.
Next, I want to thank all the present and past members in the Biomedical Ultrasonics
and Electromagnetics Lab. I thank Xiaofeng Zhao for his help on finite different simulations,
and for all of his help since the first day I joined the lab. I especially want to thank him for
being the photographer for my wedding, where he took several gorgeous photos of Qianwei
and myself. I thank Pedro Nariyoshi for his help on B-mode ultrasound simulations and for
answering a lot of my coding questions. I especially want to thank him for his humorous and
positive attitude, which makes the research life easier. I would also like to thank Yi Zhu,
Peter Beard, Drew Murray, Jim Kelly, and Kenneth Stewart. I thank all of them for the help

vi

on my projects and the interesting discussions on different aspects of medical ultrasound.
Also, I want to thank my graduate student colleagues in the ECE department. I thank
Jie Li for tolerating me for a year as a roommate, teaching me what is a Green’s function, and
providing interesting discussions on computational electromagnetics, which is not quite as
interesting as acoustics. I thank Fan Bin for interesting discussions and for his introduction
to MEMS, and I would like to thank him for his help inside and outside of graduate school.
I especially want to thank the “cool guys lunch club” (Xiaofeng Zhao, Pedro Nariyoshi, and
Fan Bin) for sharing lunch time. I hope that we will have the chance to have lunch together
in the future. I thank Xiaorui Wang, Pedro Nariyoshi, Sylmarie Davila, Patrick O’Hara,
and Yousef Gtat for working as a team in teaching the ECE 331 (microcontroller) lab. I
also want to thank Tianlong Song, Runruo Chen, Jiankun Liu, Mingwu Gao, Yu Cheng,
Yan Shi, Zhe Wang, Xiaorui Wang, and Biyi Fang for sharing knowledge in their fields that
broadened my understanding of electrical engineering.
I would like to thank Siemens Healthineers in Issaquah, WA for providing me a valuable
internship opportunity during my Ph.D. study so that I could gain some practical experience
of medical ultrasound in industry. I thank Liexiang Fan for bringing me into the advanced
development team and for sharing his experiences not only in medical ultrasound but also
in career building, I thank Yassin Labyed for his great guidance and being the best industry
mentor as you can imagine, and I thank Xiaozheng Zeng for recommending this opportunity
to me and for providing great help both in work and life.
I would like to thank all of the alumni of the Nanjing University Acoustics Alumni
Association in America (NJU-AAAA), who shared useful career suggestions with me. I
would like to offer special thanks to my undergraduate roommate Tao Sun, who is now at
Brigham and Women’s Hospital, Harvard Medical School, and Tufts University (yes, that
many affiliations), for his discussions on microbubbles and therapeutic ultrasound, and for
sharing the latest academic gossip.
I also want to thank my two cats, Ramen and Rancho, for their companionship during

vii

the hardest time of my Ph.D. study.
I want to thank my family back in China for their support and care during my Ph.D.
study, especially my parents Zhimin Peng and Xiaoning Chen and my grandmom Guozhen
Li for their love and encouragement.
Lastly and most importantly, I want to thank my wife, Qianwei Jiang, who provided
unconditional love and companionship throughout my Ph.D. study. This thesis could not
be finished without your love, support, and encouragement. I’m deeply thankful that I meet
you on orientation day for graduate school (08/16/2012) in the Breslin Center, which then
made the days of the following years the best days in my life.
Yiqun Yang
January, 2018

viii

TABLE OF CONTENTS

LIST OF TABLES

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Ultrasound elastography . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Ultrasound elastography techniques . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Static excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Natural physiological excitation . . . . . . . . . . . . . . . . . . . . .
1.2.3 Transient excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3.1 Acoustic Radiation Force Impulse (ARFI) Imaging . . . . .
1.2.3.2 Shear Wave Elasticity Imaging (SWEI) . . . . . . . . . . . .
1.2.3.3 Supersonic Shear Imaging (SSI) . . . . . . . . . . . . . . . .
1.2.3.4 Spatially Modulated Ultrasound Radiation Force (SMURF)
1.2.3.5 Comb-Push Ultrasound Shear Elastography (CUSE) . . . .
1.2.4 Harmonic excitation . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.4.1 Sonoelasticity imaging . . . . . . . . . . . . . . . . . . . . .
1.2.4.2 Vibro-Acoustography (VA) . . . . . . . . . . . . . . . . . .
1.2.4.3 Harmonic Motion Imaging (HMI) . . . . . . . . . . . . . . .
1.2.4.4 Shear wave Dispersion Ultrasound Vibrometry (SDUV) . . .
1.3 Review of numerical simulation methods for shear waves induced by an
acoustic radiation force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Acoustic radiation force simulations . . . . . . . . . . . . . . . . . . .
1.3.2 Shear wave simulations . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Motivation and thesis structure . . . . . . . . . . . . . . . . . . . . . . . . .

9
10
11
12

CHAPTER 2 BACKGROUND . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Acoustic radiation force . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Equation of motion in elastic and viscoelastic soft media . . . . . . . . . . .
2.3 Green’s function solutions to the equation of motions . . . . . . . . . . . . .
2.3.1 Green’s function for Navier’s equation in elastic media . . . . . . . .
2.3.2 Green’s function to the Navier’s equation in viscoelastic media . . . .
2.4 Shear wave measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Methods for estimating the mechanical parameters from measured shear waves
2.5.1 The time-of-flight (TOF) method . . . . . . . . . . . . . . . . . . . .
2.5.2 The k-space method . . . . . . . . . . . . . . . . . . . . . . . . . . .

14
14
14
17
17
19
21
22
22
23

CHAPTER 3 GPU-BASED GREEN’S FUNCTION SIMULATIONS OF SHEAR
WAVES GENERATED BY AN APPLIED ACOUSTIC RADIATION FORCE IN ELASTIC AND VISCOELASTIC MEDIA . . . . .
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
24

ix

1
1
2
2
3
3
4
5
5
6
6
7
7
8
8
9

3.2

3.3

3.4

3.5

3.6

Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Navier’s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Green’s functions calculations . . . . . . . . . . . . . . . . . . . . .
3.2.3 Simulating shear waves with Green’s functions . . . . . . . . . . . .
Numerical Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Acoustic radiation force . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Parallel implementation . . . . . . . . . . . . . . . . . . . . . . . .
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Simulated shear wave particle displacements at discrete observation
points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Computation times for serial and parallel calculations . . . . . . . .
3.4.3 Convergence of Green’s function calculations . . . . . . . . . . . . .
3.4.4 Comparisons between simulated and measured shear waves . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1 Diffraction, attenuation, and dispersion . . . . . . . . . . . . . . . .
3.5.2 Serial vs parallel implementations . . . . . . . . . . . . . . . . . . .
3.5.3 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.4 Comparisons between measured and simulated shear waves . . . . .
3.5.5 Intended applications . . . . . . . . . . . . . . . . . . . . . . . . . .
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

25
25
26
27
29
29
30
31

.
.
.
.
.
.
.
.
.
.
.

32
34
35
41
47
47
48
49
51
53
53

CHAPTER 4 TIME-DOMAIN ESTIMATION OF SHEAR ELASTICITY AND
SHEAR VISCOSITY FOR SHEAR WAVE ELASTICITY IMAGING (SWEI) IN VISCOELASTIC SHEAR WAVE PHANTOMS . . .
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Acoustic radiation force . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Shear waves in viscoelastic media . . . . . . . . . . . . . . . . . . . .
4.3 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Shear wave measurements . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Shear wave simulations . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Phantom experiments . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.4 Time-of-flight method for shear elasticity estimation . . . . . . . . . .
4.3.5 Model-based approach for shear elasticity and shear viscosity estimation
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Time of flight parameter estimates evaluated with simulated shear
wave data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2 New model-based shear wave parameter estimation approach applied
to simulated shear wave data . . . . . . . . . . . . . . . . . . . . . .
4.4.3 TOF and model-based estimates of shear wave parameters for viscoelastic shear wave phantoms . . . . . . . . . . . . . . . . . . . . . .
4.4.4 Shear elasticity and shear viscosity imaging in viscoelastic shear
wave phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.5 Shear elasticity and shear viscosity imaging with different push beams
4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
x

55
55
57
57
58
59
59
60
62
63
63
64
64
66
71
73
76
79

4.5.1

4.6

Comparisons between the TOF method and the new model-based
approach: simulated 3D shear wave data . . . . . . . . . . . . . . .
4.5.2 Comparisons between the TOF method and the new model-based
approach: viscoelastic shear wave phantoms . . . . . . . . . . . . .
4.5.3 Shear wave parameter estimation in different phantoms . . . . . . .
4.5.4 Shear wave parameter estimation with different push beams . . . .
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

79

.
.
.
.

80
81
83
84

CHAPTER 5 SHEAR ELASTICITY AND SHEAR VISCOSITY ESTIMATION
IN LIVER . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Acoustic radiation force . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Navier’s equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3 Green’s functions for Navier’s equation . . . . . . . . . . . . . . . . .
5.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Acoustic radiation force simulations . . . . . . . . . . . . . . . . . . .
5.3.2 Calculating shear wave particle velocities induced by an acoustic
radiation force using Green’s functions . . . . . . . . . . . . . . . . .
5.3.3 Approximate model for shear wave calculation in viscoelastic media .
5.3.4 Validation of the approximate viscoelastic model . . . . . . . . . . . .
5.3.5 Shear wave simulations with the approximate model using GPU . . .
5.3.6 Tissue mechanical parameter estimation using a genetic algorithm . .
5.3.7 Validation with In Vivo human data . . . . . . . . . . . . . . . . . .
5.4 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1 Shear wave particle velocities measured in liver . . . . . . . . . . . .
5.4.2 Shear elasticity and shear viscosity estimation with the new modelbased approach at a single point in space . . . . . . . . . . . . . . . .
5.4.3 Shear elasticity and shear viscosity imaging in human liver . . . . . .
5.4.4 Shear wave parameter estimates in liver . . . . . . . . . . . . . . . .
5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5.1 Shear waves measured with different push beams . . . . . . . . . . .
5.5.2 Shear wave parameter estimates in different locations . . . . . . . . .
5.5.3 Shear wave parameter estimates for in vivo liver . . . . . . . . . . . .
5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86
86
88
88
89
90
91
91
92
92
95
97
98
100
101
101
103
105
108
112
112
113
114
115

CHAPTER 6 CONCLUSION AND FUTURE WORK . . . . . . . . . . . . . . . . 116
6.1 Summary of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

xi

LIST OF TABLES

Table 3.1: The FWHM values of the measured and simulated shear wave particle
velocities in the x direction evaluated at the peak at 3 ms, 6 ms, and 9
ms for phantom E3962-1 and E3962-2. . . . . . . . . . . . . . . . . . . . .

51

Table 3.2: The FWHM values of the measured and simulated shear wave particle
velocities in the x direction evaluated at the peak at 3 ms, 6 ms, and 9
ms for phantom E3962-3. . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

Table 4.1: Parameters for the three push beams. . . . . . . . . . . . . . . . . . . . .

59

Table 5.1: The average and the standard deviation of the estimated shear elasticity values provided by the new model-based approach for 20 different
measurements. The units are kPa for each entry in the table. . . . . . . . 108
Table 5.2: The average and the standard deviations of the shear elasticity values
estimated with the TOF method for 20 different measurements. The
units are kPa for each entry in the table. . . . . . . . . . . . . . . . . . . 110
Table 5.3: The average and the standard deviations of the shear viscosity values
estimated with the new model-based approach for 20 different measurements. The units are Pa·s for each entry in the table. . . . . . . . . . . . 111

xii

LIST OF FIGURES

Figure 3.1: The locations of the observation points relative to a push beam focused
at 25 mm in the x-z plane with y = 0. . . . . . . . . . . . . . . . . . . .

32

Figure 3.2: Simulated shear wave particle displacements computed with elastic and
viscoelastic models for an applied 3D acoustic radiation force push evaluated at the 8 observation points indicated in Figure 3.1. The solid
lines describe the normalized shear wave particle displacements calculated for the elastic model with µ = 2.25 kPa and ηs = 0 Pa·s, and the
dashed lines describe the normalized shear wave particle displacements
calculated for the viscoelastic model with µ = 2.25 kPa and ηs = 1 Pa·s.

33

Figure 3.3: Computation times for serial and parallel implementations of time domain shear wave calculations with Green’s functions. The solid line with
circle markers indicates the computation times for parallel simulations
of shear waves with the elastic model, the dashed-dotted line with triangle markers indicates the computation times for parallel simulations of
shear waves with the viscoelastic model, the dashed line with diamond
markers indicates the computation times for serial simulations of shear
waves with the elastic model, and the dotted line with square markers
indicates the computation times for serial simulations of shear waves
with the viscoelastic model. . . . . . . . . . . . . . . . . . . . . . . . . .

36

Figure 3.4: 2D cross-sections of 3D ultrasound beams simulated in FOCUS. Figures
(3.4a), (3.4b), (3.4c), and (3.4d) show the push beams with f/1.5, f/2,
f/3, and f/∞, respectively. In (3.4a)-(3.4d), each beam is normalized
by the maximum intensity value. . . . . . . . . . . . . . . . . . . . . . .

37

Figure 3.5: Percent differences calculated for shear wave particle displacements simulated with the elastic model. The shear elasticity and shear viscosity
values in these simulations are µ = 2.25 kPa and ηs = 0 Pa·s, where the
push beams are calculated with different spatial samplings. The solid
line with circle markers indicates the percent difference for the shear
wave particle displacements induced by an unfocused push beam with
f/∞, the dotted line with square markers indicates the percent difference for the shear wave particle displacements induced by a focused
push beam with f/3, the dash-dotted line with triangle markers indicates the percent difference for the shear wave particle displacements
induced by a focused push beam with f/2, and the dashed line with
diamond markers indicates the percent difference for the shear wave
particle displacements induced by a focused push beam with f/1.5. . . . .

38

xiii

Figure 3.6: Percent differences calculated for shear wave particle displacements simulated with the elastic model. The shear elasticity and shear viscosity
values in these simulations are µ = 2.25 kPa and ηs = 1 Pa·s, where
the push beams are calculated with different spatial samplings. The
descriptions of the markers and line types are the same as those given
in the caption for Figure 3.5. . . . . . . . . . . . . . . . . . . . . . . . . .

40

Figure 3.7: Measured shear wave particle velocities in phantom E3962-1 and the
corresponding simulated shear wave particle velocities that minimize
Eq. 3.22 at three different times. Figs. 3.7(a), 3.7(b), and 3.7(c) show
the shear wave particle velocities measured at 3 ms, 6 ms, and 9 ms,
respectively. Figs. 3.7(d), 3.7(e), and 3.7(f) show the corresponding
simulated shear wave particle velocities computed with µ = 0.3 kPa
and ηs = 0.1 Pa·s at 3 ms, 6 ms, and 9 ms, respectively. . . . . . . . . .

43

Figure 3.8: Measured shear wave particle velocities in phantom E3962-2 and the
corresponding simulated shear wave particle velocities that minimize
Eq. 3.22 at three different times. Figs. 3.8(a), 3.8(b), and 3.8(c) show
the shear wave particle velocities measured at 3 ms, 6 ms, and 9 ms,
respectively. Figs. 3.8(d), 3.8(e), and 3.8(f) show the corresponding
simulated shear wave particle velocities with µ = 2.3 kPa and ηs = 0.1
Pa·s at 3 ms, 6 ms, and 9 ms, respectively. . . . . . . . . . . . . . . . . .

44

Figure 3.9: Measured shear wave particle velocities in phantom E3962-3 and the
corresponding simulated shear wave particle velocities that minimize
Eq. 3.22 evaluated at three different times. Figs. 3.9(a), 3.9(b), and
3.9(c) show the shear wave particle velocities measured at 2 ms, 4 ms,
and 6 ms, respectively. Figs. 3.9(d), 3.9(e), and 3.9(f) show the corresponding simulated shear wave particle velocities computed with µ =
4.9 kPa and ηs = 1.5 Pa·s at 2 ms, and computed with µ = 4.5 kPa and
ηs = 1.5 Pa·s at 3 ms and 4 ms, respectively. . . . . . . . . . . . . . . . .

46

Figure 4.1: 2D cross-sections of the three push beams simulated in FOCUS. Figures (4.1a), (4.1b), and (4.1c) show push beams with focal depths of 33
mm, 45 mm, and 70 mm, respectively. . . . . . . . . . . . . . . . . . . .

60

Figure 4.2: The measured shear wave particle velocity induced by the push beam
with a focal depth of 33 mm in phantom E2348-1. Figures (4.2a), (4.2b),
and (4.2c) show the measured shear wave particle velocity at 2 ms, 6
ms, and 10 ms, respectively. . . . . . . . . . . . . . . . . . . . . . . . . .

62

xiv

Figure 4.3: Shear elasticity images generated by the TOF method using noisy 3D
shear wave particle velocity data obtained from elastic and viscoelastic
simulation models. Fig. 4.3a) shows the estimated shear elasticity for
the homogeneous elastic simulation model with µ = 4 kPa and ηs =
0 Pa·s, and Fig. 4.3b) shows the estimated shear elasticity for the homogeneous viscoelastic simulation model with µ = 4 kPa and ηs = 2
Pa·s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

Figure 4.4: Objective functions for the shear wave particle velocities calculated at
(x, z) = (10 mm, 30 mm) with the elastic and viscoelastic simulation
models. Fig. 4.4a shows the objective function evaluated for the elastic
simulation model, and Fig. 4.4b shows the objective function evaluated
for the viscoelastic simulation model. The “o” markers indicate the
estimates obtained with the new model-based approach, where µ = 4.00
kPa and ηs = 0 in Fig. 4.4a and µ = 4.00 kPa and ηs = 1.90 Pa·s in
Fig. 4.4b. The “x” markers indicate the estimated values obtained with
the TOF method, where µ = 4.11 kPa and ηs = 0 in Fig. 4.4a and µ =
5.69 kPa and ηs = 0 in Fig. 4.4b. . . . . . . . . . . . . . . . . . . . . . .

68

Figure 4.5: Shear elasticity and shear viscosity images generated with parameters
obtained from the new model-based approach applied to simulated noisy
shear wave particle velocities calculated with the elastic and viscoelastic
simulation models. Figs. 4.5(a) and 4.5(b) show the shear elasticity and
shear viscosity images, respectively, for the elastic simulation model
with µ = 4 kPa and ηs = 0 Pa·s, respectively. Figs. 4.5(c) and 4.5(d)
show the shear elasticity and shear viscosity images, respectively, for a
viscoelastic simulation model with µ = 4 kPa and ηs = 2 Pa·s. . . . . . .

69

Figure 4.6: The objective function in Eq. 4.9 evaluated for different values of the
shear elasticity and shear viscosity for shear wave particle velocities
measured at (x, z) = (10 mm, 30 mm) in phantom E2348-1. The “o”
marker indicates the value obtained with the new model-based estimation approach (µ = 3.60 kPa and ηs = 0.80 Pa·s), and the “x” marker
indicates the value produced by the TOF method (µ = 6.70 kPa and ηs
= 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

xv

Figure 4.7: Measured and simulated shear wave particle velocities using shear wave
parameter values estimated with the time-of-flight (TOF) method and
the new model-based estimation approach. In Figs. 4.7a and 4.7b, the
blue dashed line indicates the measured shear wave particle velocity
induced by the push beam with a focal peak at 33 mm in phantom
E2348-1, the black dotted line in Fig. 4.7(a) is the the simulated shear
wave particle velocity with the shear wave parameters estimated with
the TOF method (µ = 6.70 kPa and ηs = 0), and the red solid line
in Fig. 4.7(b) is the the simulated shear wave particle velocity with the
shear wave parameters (µ = 3.60 kPa and ηs = 0.80 Pa·s) obtained from
the new model-based estimation approach. . . . . . . . . . . . . . . . .

72

Figure 4.8: Shear wave parameter images generated by the TOF method and the
new model-based approach for phantoms E2348-1, E2348-2, and E23483 where the shear wave is induced by a push beam focused at 45 mm.
Figs. 4.8(a)-(c) show the shear elasticity image provided by the TOF
method, the shear elasticity image provided by the new model-based approach, and the shear viscosity image provided by the new model-based
approach evaluated for phantom E2348-1, respectively. Figs. 4.8(d)-(f)
show the shear elasticity image provided by the TOF method, the shear
elasticity image provided by the new model-based approach, and the
shear viscosity image provided by the new model-based approach evaluated for phantom E2348-2, respectively. Figs. 4.8(g)-(i) show the shear
elasticity image provided by the TOF method, the shear elasticity image
provided by the new model-based approach, and the shear viscosity image provided by the new model-based approach evaluated for phantom
E2348-3, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

Figure 4.9: Shear wave parameter images generated by the TOF method and the
new model-based approach for phantom E2348-2. Figs. 4.9(a), 4.9(e),
and 4.9(i) show the shear elasticity images generated by the TOF method
with push beams focused at 33 mm, 45 mm, and 70 mm, respectively.
Figs. 4.9(b), 4.9(f), and 4.9(j) show the shear elasticity images generated
by the new model-based approach with push beams focused at 33 mm,
45 mm, and 70 mm, respectively. Figs. 4.9(c), 4.9(g), and 4.9(k) show
the shear viscosity images generated by the new model-based approach
with push beams focused at 33 mm, 45 mm, and 70 mm, respectively.
Figs. 4.9(d), 4.9(h), and 4.9(l) show the images of the difference between
the measured shear wave particle velocity and the simulated shear wave
particle velocity for push beams focused at 33 mm, 45 mm, and 70 mm,
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

xvi

Figure 4.10: Comparison of the shear elasticity values estimated in three CIRS viscoelastic shear wave phantoms with the new model-based approach and
the time-of-flight method. The three groups of bars represent the results for three different phantoms. For each group of bars, the first three
bars indicate the result obtained with the new model-based approach,
and the last three bars represent the result given by the time-of-flight
method. The error bars represent the standard deviation of the estimated elasticity calculated from ten measurements. . . . . . . . . . . . .

82

Figure 4.11: Estimated viscosity values obtained with the model-based approach in
three CIRS viscoelastic shear wave phantoms. Three different groups
of bars describe the results for three different phantoms. The three
bars in each group indicate the results generated by pushes at different
depths. The error bars indicate the standard deviation of the estimated
elasticity from ten measurements. . . . . . . . . . . . . . . . . . . . . . .

83

Figure 5.1: The results of the individual calculation steps that produce an approximate shear wave particle displacement. Fig. 5.1a shows the push beam,
and the red dot indicates the observation point, which is located at
(x, z) = (10 mm, 22.8 mm). Fig. 5.1b shows the shear wave particle
displacement simulated with the elastic model Eq. 5.22. Fig. 5.1c shows
the attenuation term calculated with Eq. 5.18. Fig. 5.1d shows the approximate shear wave particle displacement calculated by convolving
the results shown in Fig. 5.1b and Fig. 5.1c. . . . . . . . . . . . . . . . .

95

Figure 5.2: Comparison between the shear wave particle displacements simulated
with the 3D viscoelastic model and approximate model. The blue lines
indicate the shear wave particle displacements calculated with the viscoelastic model, and the red dashed lines represent shear wave particle
displacements computed with the approximate model. . . . . . . . . . . .

97

Figure 5.3: A 2D image showing the percent difference between shear wave particle
displacements calculated with the approximate model in Eq. 5.21 and
the viscoelastic model in Eqs. 5.14. The black contour represents a value
of 5% for the percent difference. . . . . . . . . . . . . . . . . . . . . . . .

98

Figure 5.4: Two push beams that induce shear waves in human liver. Fig. 5.4a
shows the push beam focused at 35 mm, and Fig. 5.4b shows the push
beam focused at 50 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Figure 5.5: The shear wave particle velocities measured in human liver. Figs. 5.5(a)5.5(c) show the shear wave particle velocities induced by the push beam
focused at 35 mm at 2 ms, 4 ms, and 6 ms, respectively, and Figs. 5.5(d)5.5(f) show the shear wave particle velocities induced by the push beam
focused at 50 mm at 2 ms, 4 ms, and 6 ms, respectively. . . . . . . . . . 102
xvii

Figure 5.6: Comparisons between measured shear wave particle velocities and shear
wave particle velocities simulated with the µ and ηs values estimated
using the new model-based approach at different locations in space. The
blue solid lines indicate the measured shear wave particle velocities,
and the red dashed lines represent the simulated shear wave particle
velocities. Fig. 5.6(a) shows the measured and simulated shear wave
particle velocities at (x, z) = (5 mm, 30 mm), Fig. 5.6(b) shows the
measured and simulated shear wave particle velocities at (x, z) = (6
mm, 35 mm), Fig. 5.6(c) shows the measured and simulated shear wave
particle velocities at (x, z) = (3 mm, 30 mm), and Fig. 5.6(d) shows the
measured and simulated shear wave particle velocities at (x, z) = (15
mm, 30 mm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Figure 5.7: Shear wave parameter images for human liver produced by a push beam
focused at 35 mm. Fig. 5.7(a) shows the shear elasticity image generated
by the TOF method. Fig. 5.7(b) shows the shear elasticity image generated by the new model-based approach. Fig. 5.7(c) shows the shear
viscosity image generated by the new model-based approach. Fig. 5.7(d)
shows the percent difference between the measured shear waves and the
simulated shear waves with the µ and ηs values estimated with the new
model-based approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
Figure 5.8: Shear wave parameter images for human liver produced by a push beam
focused at 50 mm. Fig. 5.8(a) shows the shear elasticity image generated
by the TOF method. Fig. 5.8(b) shows the shear elasticity image generated by the new model-based approach. Fig. 5.8(c) shows the shear
viscosity image generated by the new model-based approach. Fig. 5.8(d)
shows the percent difference between the measured shear waves and the
simulated shear waves with the µ and ηs values estimated with the new
model-based approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

xviii

CHAPTER 1
INTRODUCTION

1.1

Ultrasound elastography

For many years, physicians have used the stiffness of biological tissues to identify the
pathological conditions of organs to diagnose diseases. As early as 400 B.C., Hippocrates
discovered that the normal abdominals are “soft, free of pain, and yield when pressed with the
finger, are more chronic,” while abdominal swellings are “painful, hard, and large, indicate
danger of speedy death” [1]. Clinically, palpation is used to qualitatively detect the stiffness
of the tissue for early cancer screening in breast [2] or liver [3]. However, although palpation
is a convenient and straightforward diagnostic approach, it has several significant limitations.
First, it is only useful in examining shallower tissue, and the effectiveness of this technique
drops dramatically when the tissue is deep inside the body. Also, manual palpation is a
qualitative technique, so no quantitative information for tissue can be obtained. Thus, the
diagnostic accuracy relies on the experience of the physician. All of the above limitations
suggest that better diagnostic techniques are needed to assess the stiffness of the tissue.
Ultrasound elastography, which has been actively developed for the past two decades, is
an enhanced palpation technique that overcomes these limitations. Ultrasound elastography
performs “remote palpation,” which noninvasively detects the mechanical properties of the
tissue deep inside the body using ultrasound. The procedure for ultrasound elastography is
generally described as [4–10]:
1. Induce a tissue deformation with an external force. The external force can be a
static compression force, a natural physiological force such as a pulsating blood vessel, or an
acoustic radiation force (ARF) that is generated by a high-intensity ultrasound beam.
2. Measure the tissue deformation using ultrasound. Typically, the measurement is
performed by comparing the ultrasound pulse-echo signal before and after the application

1

of the external force. Correlation-typed algorithms are applied to detect the movement of
local ultrasound speckle [11] and then measure the local movement of tissue.
3. Estimate the tissue mechanical properties using the models that relate the measured
tissue deformation to the mechanical properties.
In this introductory chapter, multiple ultrasound elastography techniques are reviewed
based on different types of excitations. In addition, methods that numerically simulate
ultrasound elastography are reviewed. Finally, the motivation for this work is described and
the outline of this thesis is given.

1.2

Ultrasound elastography techniques

Ultrasound elastography techniques can be classified based on the type of the external
force that induces the tissue deformation. In general, four types of excitation sources are
used to induce the deformation in the tissue: static excitation, natural physiological excitation, transient excitation, and harmonic excitation. In this section, multiple ultrasound
elastography techniques based on these four types of excitation schemes are reviewed.

1.2.1

Static excitation

With static excitation, tissue is deformed by a static external force, which is usually a manual
compression applied by the sonographer. A representative example is elastography, which
was introduced by Ophir et al. in 1991 [6,12]. In this technique, ultrasound measurement of
tissue is first performed at the equilibrium state with no external force applied. Then, the
tissue is compressed manually, and ultrasound measurements of the tissue are taken when the
tissue reaches a static compressed status. Comparisons are made between the ultrasound
pulse-echo signals measured before and after the tissue is compressed, and the images of
the local tissue strain, which are also called “elastograms,” are generated. Typically, lower
strains correspond to harder tissue, and higher strains correspond to softer tissue. Clinically, this technique is used for assessing the hardness of muscle [13–16], detecting breast

2

tumors [13,17–20], detecting thyroid tumors [21–23], detecting prostate tumors [24–26], characterizing thermal ablation lesions after high-intensity focused ultrasound (HIFU) or radio
frequency (RF) treatments [27–29], and detecting vascular plaques with a static excitation
deployed by an intravascular balloon [30,31]. The advantage of this technique is that the excitation is relatively easy to apply, but the main disadvantage is that only qualitative images
that differentiate between tissues with different stiffnesses are generated but no quantitative
information describing the tissue stiffness is provided.

1.2.2

Natural physiological excitation

Natural physiological motion, such as a pulsating blood vessel, heartbeat, or muscle tension
can also induce tissue deformation, and ultrasound elastography techniques that exploit this
behavior have been developed to assess tissue stiffness. In these techniques, the method
for measuring tissue deformation is similar to ultrasound elastography with a static excitation. Ultrasound measurements are performed at different times, while the physiological
excitation occurs intermittently, such as the periodic pulsing of a blood vessel. Comparisons
are performed between the ultrasound pulse-echo signals measured before and after tissue
deformation, and elastograms are generated based on the amplitude of the strain. Clinically,
these methods assess the stiffness of vessels and detect vascular plaques [32–35], characterize
the stiffness of muscle [36, 37], and detect the stiffness of myocardial tissue [38–40]. The advantage of this technique is that no manual or external excitation is needed to induce tissue
deformation, but the obvious drawback is that it is not possible to control the excitation. In
addition, this technique cannot provide quantitative information about tissue stiffness.

1.2.3

Transient excitation

In transient ultrasound elastography, tissue deformation is induced using an acoustic radiation force (ARF) for a short period of time (typically several µs to several hundreds
of µs) [41, 42]. The ARF is generated by a high-intensity ultrasound beam that is often
3

called the “push beam”. Compared to ultrasound elastography techniques based on static
compression or natural physiological motion, ARF-based ultrasound elastography techniques
have several advantages including a controllable excitation, so different tissue deformation
responses can be induced by changing the type of ARF excitation. Furthermore, ARF can
be applied to tissue deep inside the body, which expands its applications beyond superficial targets. Moreover, shear waves can be induced by an ARF push beam in tissue, and
quantitative mechanical properties of tissue can be obtained by measuring the shear wave
propagation speed. A limitation of this technique is that useful information is only obtained
around the focus of the push beam, but this limitation is overcome by applying multiple push
beams focused at different depths. Multiple variations of transient ultrasound elastography
are reviewed in sections 1.2.3.1-1.2.3.5.

1.2.3.1

Acoustic Radiation Force Impulse (ARFI) Imaging

Acoustic Radiation Force Impulse (ARFI) Imaging technique was introduced by Nightingale
et al. [41] in 2000. In ARFI, an ultrasound transducer array generates a focused ARF push
beam, and the push beam produces tissue deformation. Before and after the push beam is
applied to the tissue, multiple frames of radio-frequency (RF) or in-phase/quadrature (IQ)
data at the location of the push beam are obtained with pulse-echo ultrasound sequences.
Then, the particle displacement is calculated from the IQ data using one-dimensional (1D)
autocorrelations between the data measured before and after the push is applied [11, 43].
Then, the push beams and the corresponding ultrasound measurements are scanned through
the whole two-dimensional (2D) field-of-view (FOV), and qualitative images for particle
displacement are generated. Typically, a large displacement suggests that the tissue is soft,
while a small displacement suggests that the tissue is hard. In clinical practice, ARFI imaging
has been used to examine tissue stiffness and to detect tumors in liver [44,45], thyroid [46,47],
breast [48,49], and prostate [50,51]. ARFI has also been used to assess elasticity properties in
vascular tissue [52–54], detecting plaque in vascular tissue [55, 56], and characterize stiffness

4

in cardiac tissues [57, 58]. Moreover, ARFI has been used to guide thermal ablation therapy
by monitoring the change in stiffness of the tissue during treatments [59]. ARFI imaging
only provides qualitative information that describes the relative stiffness, but no quantitative
information describing the tissue elasticity is given.

1.2.3.2

Shear Wave Elasticity Imaging (SWEI)

The shear wave elasticity imaging (SWEI) technique was introduced by Sarvazyan et al. in
1998 [42]. In SWEI, one or more focused or unfocused ultrasound beams are applied for
a short period of time to generate tissue deformation. In turn, this deformation produces
shear waves. To measure the shear wave particle displacements, multiple frames of RF or
IQ data within a 2D FOV are obtained with pulse-echo ultrasound sequences. Then, shear
wave particle displacements are calculated from the RF or IQ data using 1D autocorrelations between consecutive frames [11, 42, 60]. Next, the shear wave speed is estimated from
the measured particle displacements, and quantitative shear elasticity parameters are then
obtained from the estimated shear wave speed. In clinical practice, SWEI has been used
to quantitatively measure the stiffness of liver tissue and determine the stage of liver fibrosis [44,61–66]. Other clinical applications for SWEI include quantitative assessment of tissue
mechanical properties in thyroid [67, 68], breast [69, 70], kidney [71], and prostate [72, 73].

1.2.3.3

Supersonic Shear Imaging (SSI)

The supersonic shear imaging (SSI) technique was introduced by Bercoff et al. in 2004 [74].
SSI is similar to SWEI, but the shear wave generation and detection approaches are different.
In SSI, several push beams focused at different depths are sequentially applied to the tissue.
The focus of the push beams moves downward faster than the propagation speed of the
shear waves, hence the term “supersonic.” After the sequence of push beams is applied,
shear waves are then generated that propagate away from the push location with a coneshaped wavefront. After the shear waves are generated, multiple frames of RF or IQ data
5

within a 2D FOV are measured using an ultrafast compound imaging technique [75, 76],
which measures the ultrasound signals at a very high frame rate (up to 6000 frames per
second). Then, shear wave particle displacements are calculated using 1D autocorrelations
between consecutive frames [11]. In clinical practice, SSI has been used to quantitatively
assess the shear elasticity of liver [77–80], breast [81, 82], cornea [78], muscle [83–85], and to
guide thermal ablation [74].

1.2.3.4

Spatially Modulated Ultrasound Radiation Force (SMURF)

The spatially modulated ultrasound radiation force (SMURF) technique was developed by
McAleavey et al. in 2007 [86]. As apposed to SWEI and SSI, which excite the tissue in a
fixed location and then measure the signal at different lateral positions, the SMURF technique induces shear waves at different lateral positions by spatially moderating the acoustic
radiation force transversely. SMURF then measures the induced shear waves at the same
position. By measuring the arrival time difference of the shear waves generated by these
different pushes located at known positions, the shear wave speed is estimated and the elasticity parameters are obtained. The SMURF technique suffers from fewer problems with
correlated noise compared with other transient ultrasound elastography techniques because
shear waves are measured at the same location, so the noise and jitter are reduced [87].

1.2.3.5

Comb-Push Ultrasound Shear Elastography (CUSE)

The comb-push ultrasound shear elastography (CUSE) technique was developed by Song
et al. in 2012 [88]. In CUSE, the push beam is a comb-shaped acoustic radiation force
consisting of several unfocused [89] or focused [90] ultrasound beams distributed at equal
intervals within the tissue. Each push beam is generated by an individual sub-aperture of
a transducer array. The array that generates the push beam is also used as the imaging
transducer that measures the shear waves along A-lines outside of the push beam. CUSE
provides an estimate of the shear wave speed for the entire field-of-view (FOV) simultaneously
6

within a short period of time (less than 35 ms). Thus, this approach provides images of tissue
mechanical properties at a relatively high frame rate (up to approximately 30 fps) [90]. In
clinical practice, CUSE has been applied to quantitative measurements of elasticity and to
detect tumors in breast [91] and liver [92].

1.2.4

Harmonic excitation

In ultrasound elastography with harmonic excitations, a continuous wave (CW) pressure
field with frequency between 10 Hz to 1 kHz is applied to the tissue [10, 93]. The CW
excitation can be applied either with a mechanical vibrator [93] or with a special distribution
of ARF pushes [94, 95]. Compared with some other ultrasound elastography techniques, the
advantage of harmonic ultrasound elastography is that quantitative information about tissue
mechanical properties in the frequency domain can be obtained. However, a disadvantage is
that the acquisition time to obtain the mechanical properties for the entire 2D FOV is much
longer than for other ultrasound elastography techniques, so generating a full 2D image using
this technique can be challenging. Multiple variations of harmonic ultrasound elastography
are reviewed in sections 1.2.4.1-1.2.4.4.

1.2.4.1

Sonoelasticity imaging

Sonoelasticity imaging was developed by Lerner et al. in 1990 [93, 96]. In sonoelasticity
imaging, a mechanical vibrator applies an external vibration to generate oscillations at low
frequencies (10 Hz - 1 kHz) in tissue. The vibratory tissue motion is detected using Doppler
ultrasound, and the tissue mechanical properties are obtained using model-based approaches
by analyzing tissue motion at different frequencies. The drawback of this technique is that a
vibrating source is required to induce tissue deformation, so applications are depth-limited.
In clinical practice, sonoelasticity imaging assesses tissue mechanical properties in muscle
[97] and detects tumors in the liver [98, 99], kidney [98, 100], breast [98, 101], and prostate
[102, 103].
7

1.2.4.2

Vibro-Acoustography (VA)

Vibro-acoustography (VA) was developed by Fatemi et al. in 1998 [94, 104]. In VA, two
ARF push beams are simultaneously focused at the same spot with slightly different driving
frequencies. The frequency difference, or beat frequency, is on the order of 1 kHz. For an
ARF push applied in this manner, a sinusoidal vibration at the beat frequency is induced
at the focal location, and an acoustic emission is generated by the vibration. The acoustic
emission is measured using a hydrophone, and the tissue mechanical properties are obtained
by analyzing the amplitude and phase of the measured acoustic emission. The drawbacks of
this technique are that VA requires more than one transducer to transmit the excitation and
that an additional hydrophone is required to measure the response. However, the former
problem can be mitigated by using a single amplitude-modulated (AM) ARF push beam
that produces a similar vibration [95]. In addition, the long acquisition time also limits
potential in clinical applications that generate a full 2D image. In clinical practice, VA has
been used for tissue mechanical imaging in the liver [105], breast [106,107], thyroid [108,109],
prostate [110, 111], and in arteries [112].

1.2.4.3

Harmonic Motion Imaging (HMI)

Harmonic motion imaging (HMI) was developed by Konofagou and Hynynen in 2003 [113].
HMI uses the same technique as VA to generate tissue vibration, which applies two ARF
push beams focused at the same spot with slightly different driving frequencies or a single
amplitude-modulated (AM) ARF push beam. However, the technique for tissue vibration
measurement in HMI is different from that of VA. In HMI, another imaging ultrasound
transducer is used to measure the vibration instead of the hydrophone that is used in VA. The
tissue vibration is measured using pulse-echo ultrasound combined with a cross-correlation
speckle tracking algorithm [11], which is similar to the approach used in most transient
ultrasound elastography techniques. Tissue mechanical properties are obtained by analyzing
the measured tissue vibration patterns. Although HMI does not require a hydrophone, the
8

long acquisition time associated with HMI is still a disadvantage. In clinical practice, HMI
has been used to monitor thermal ablation with high intensity focused ultrasound, where the
combined approach for imaging and therapy is called harmonic motion imaging for focused
ultrasound (HMIFU) [114].

1.2.4.4

Shear wave Dispersion Ultrasound Vibrometry (SDUV)

Shear wave dispersion ultrasound vibrometry (SDUV) was developed by Chen et al. in
2004 [115, 116]. SDUV uses the same excitation method as VA to generate harmonic tissue
motion at a specific frequency. However, rather than measuring the acoustic emission from
the tissue vibration using a hydrophone, the shear wave induced by the harmonic tissue
motion is measured instead. The measurement procedure is similar to that of SWEI, where
the shear wave is measured using pulse-echo ultrasound combined with a cross-correlation
speckle tracking algorithm [11]. The shear wave propagation speed is measured by tracking the phase shift of the time-harmonic shear wave at different locations. With different
driving frequencies, the shear wave propagation speed at multiple frequencies is obtained,
and model-based [116] or model-free [117] approaches are used to estimate the tissue mechanical properties. This technique, which is performed with a single commercial ultrasound
transducer array, provides quantative measurements of shear wave dispersion for soft tissue.
In clinical applications, SDUV has been used to estimate the elasticity parameters in muscle [118], liver [116,119], kidney [120,121], prostate [122], heart [123], and blood vessels [124].

1.3

Review of numerical simulation methods for shear waves induced by an acoustic radiation force

Numerical simulations provide information that is valuable for validating and improving
ultrasound elastography techniques. In this section, a review of numerical simulation methods for transient ultrasound elastography is provided. These simulations consist of two main
parts. The first part calculates the acoustic radiation force, and the second part computes

9

the shear waves induced by the acoustic radiation force.

1.3.1

Acoustic radiation force simulations

Several ultrasound simulation toolboxes such as FOCUS, the “Fast Object-oriented
C++ Ultrasound Simulator” (http://www.egr.msu.edu/∼fultras-web/) [125], Field II
(http://field-ii.dk/) [126], Ultrasim (http://heim.ifi.uio.no/∼ultrasim/) [127], and DREAM
(http://www.signal.uu.se/Toolbox/dream/) [128] are able to simulate a steady-state focused
ultrasound beam. From the output of these, the approximate acoustic radiation force is calculated. Among the above simulation toolboxes, Field II and FOCUS are most widely used
by the medical ultrasound research community. In Field II, the ultrasound pressure field
is calculated by evaluating the impulse response of the transducer [129]. Because Field II
uses a far-field approximation, the calculation for a large transducer element is performed
by subdividing the transducer element and superposing the contribution from each small
sub-element. In FOCUS, the fast nearfield method (FNM) computes an input pressure
plane (usually λ/4 from the front face of the transducer, where λ is the wavelength of the
compressional wave) [130, 131], and then the angular spectrum approach (ASA) determines
the pressure field throughout the remaining 3D volume [132, 133]. Because the FNM does
not utilize the far-field approximation, no subdividing of elements is required. The FNM
achieves better accuracy in much less time than the impulse response method [134–136].
Furthermore, the calculation time in FOCUS for calculating ultrasound pressure fields in
three-dimensions (3D) is significantly smaller than that for impulse response methods [135],
since the combination of the FNM and ASA approaches are especially suitable for volumetric calculations. In this thesis, all of the acoustic radiation force push beams are calculated
using the FNM combined with the ASA in FOCUS.

10

1.3.2

Shear wave simulations

At present, most of the numerical simulations of shear waves generated by an applied acoustic radiation force employ the finite element (FEM) method. In 2000, Nightingale et al.
developed a finite element model which simulates the tissue response to an acoustic radiation force in breast tumors [137]. In their model, the ARF push beam is calculated using
Field II, and the shear wave induced by the push beam is calculated using LS-DYNA3D
(Livermore Software Technology Corporation, Livermore, CA). The FEM models are also
used to simulate the tissue response in ARFI imaging [138], quantifying the hepatic shear
modulus [66], demonstrating the feasibility of ARFI imaging in abdomen examinations [139],
modeling the shear wave propagation induced by an acoustic radiation force in a viscoelastic
medium [140], and other applications [139, 141, 142]. These FEM models successfully simulate shear waves induced by an acoustic radiation force, which are useful for developing
various ultrasound elastography techniques. The main drawbacks of FEM simulations are
that FEM models are usually complicated and the computation times are typically very long.
Other methods have been developed to simulate shear waves induced by acoustic radiation force push beams. The finite difference (FD) method, which is often easier to implement
than the FEM method, is also useful for simulating shear waves [143–145]. The finite difference method is sometimes faster than the FEM method, but the ability of finite difference
calculations to manage complicated boundaries and peculiar shapes is usually not as good
as FEM. The pseudospectral method is also used to simulate shear waves induced by an
acoustic radiation force [146]. In the pseudospectral method, much coarser spatial samplings
are possible for homogeneous problems relative to FD and FEM methods, which can reduce
computation time and memory usage.
The Green’s function approach is another method for simulating shear waves [147, 148].
The Green’s function approach, which is applicable to shear wave modeling in elastic [147]
and viscoelastic media [148], is appropriate when calculations at only a few observation
points in space are needed, as opposed to other numerical methods (FEM, FD, and pseu-

11

dospectral) that require the calculation of an entire 3D field. Thus, shear wave calculations
with Green’s functions can be much faster than shear wave calculations with other numerical methods in certain circumstances. Furthermore, Green’s function calculations are very
easily parallelized, since calculations for different observation points are independent of one
another. The main disadvantage is that the Green’s function approach can only simulate
shear wave propagation in homogeneous media.

1.4

Motivation and thesis structure

Although significant progress has been made in ultrasound elastography over the past
two decades, a substantial amount of work still needed to improve the image quality for
ultrasound elastography. Currently, one important difficulty in ultrasound elastography is
caused by the viscoelasticity of soft tissue. Most models that are currently used in ultrasound
elastography assume that tissue is lossless, even though multiple studies indicate that tissue
is always lossy (viscoelastic) [79, 84, 116, 121, 149–151]. Thus, estimation bias is introduced
because of this approximation [79, 116]. Quantitative ultrasound methods such as VA [94,
104–112] or SDUV [115–124] are able to characterize the viscoelasticity of tissue. However,
these methods are limited by long acquisition times and specific hardware requirements. A
robust technique that provides accurate quantitative images of viscoelastic tissue properties
with a conventional ultrasound system is a critical unmet need in medical ultrasound.
To address this problem, this thesis introduces a new approach that generates quantitative
images of viscoelastic soft tissue. Rather than analyzing the shear wave in the frequencydomain, as in time-harmonic ultrasound elastography, this new approach estimates tissue
mechanical parameters using shear waves that are measured in the time-domain using exactly
the same configuration that is employed in traditional transient ultrasound elastography (i.e.
SWEI). Thus, this new approach utilizes the same framework that is employed in current
commercial ultrasound elastography systems. In this new approach, shear waves are first
simulated using a Green’s function-based method that evaluates the Voigt model on a GPU,

12

and then the shear elasticity and shear viscosity values are estimated using an optimizationbased approach by comparing measured shear wave particle velocities with simulated shear
wave particle velocities in the time-domain. This comparison is then performed on a pointby-point basis to generate images. The new approach is evaluated with shear wave data
obtained from numerical simulations, viscoelastic shear wave phantoms, and in vivo human
livers, and the results are compared with the results of other methods for estimating shear
wave parameters.
The structure of this thesis is described as follows. In Chapter 2, background information is introduced on acoustic radiation force, elastodynamics, shear wave measurement
techniques, and shear wave elasticity estimation methods. In Chapter 3, a GPU-based
numerical simulation approach for the calculation of shear waves induced by an acoustic
radiation force in elastic and viscoelastic soft tissue models is described. Examples of shear
waves simulated using this approach are shown, and the numerical simulation methods are
validated by comparing the simulated shear wave particle velocities with shear wave particle velocities measured in several different viscoelastic shear wave phantoms. In Chapter
4, an model-based approach that estimates shear elasticity and shear viscosity values by
minimizing the difference between measured and simulated shear wave particle velocities is
introduced. The model-based estimation approach is validated with results obtained from
numerical simulations, and the model-based estimation approach is also used to estimate
shear wave parameters from measured shear wave data collected in viscoelastic shear wave
phantoms. Results suggest that the new model-based estimation approach yields more accurate parameter values than the conventional TOF method. In Chapter 5, the new approach
is accelerated with an approximation, and then the new approach generates estimates of
shear wave parameters with measured data obtained from in vivo human livers. In Chapter
6, the summary and potential future directions are presented.

13

CHAPTER 2
BACKGROUND

2.1

Acoustic radiation force

In clinical ultrasound elastography, shear waves in soft tissue are induced by an applied
acoustic radiation force, which is generated from a focused ultrasound beam. The acoustic
radiation force is modeled as [41]
f (r, t) =

2αI(r, t)
,
cp

(2.1)

where f (r, t) is the body force, r denotes a spatial point in Cartesian coordinates, t denotes
the time, α is the attenuation coefficient of the compressional wave, cp is the compressional
wave speed, and I(r, t) is the acoustic intensity. The intensity of the ultrasound beam is
proportional to the square of the acoustic pressure p, which is defined as I = |p|2 /(ρcp ).
In transient elastography, the acoustic radiation force is applied to the tissue for a short
period of time to produce a shear wave. To account for the temporal component of the
acoustic radiation force, f (r, t) is written as f (r, t) = f (r)w(t), where f (r) is the spatial
distribution of the acoustic radiation force, and w(t) is a window function that describes
the time dependence of the acoustic radiation force. The window function is defined here as
w(t) = rect[(t − T /2)/T ], where rect(·) is the rectangle function, and T is the duration of
the push.

2.2

Equation of motion in elastic and viscoelastic soft media

The elastodynamic response induced by an applied acoustic radiation force in elastic and
viscoelastic media is described by the equation of motion, which can be found in any classical
mechanics textbook such as [152, 153]. The derivation of the equation of motion is based
on three fundamental relations: the conservation of momentum, the strain-displacement
relations, and the stress-strain relationship. These three relationships are described below.
14

Newton’s second law guarantees the conservation of momentum, which relates the stress
and the external force to the acceleration. The conservation of momentum equation is written
as [152]
∂σij
∂ 2u
+ fi = ρ 2i ,
∂xj
∂t

(2.2)

where σij is the stress tensor, fi is the external body force, ui is the displacement, i, j = 1, 2, 3
are the indices for all three dimensions, t is the time, and ρ is the density of the material.
Assuming that the strain and displacement are both small so that the behavior is linear, the
strain-displacement relationship is written as [152]
1
εij =
2



∂ui ∂uj
+
∂xj
∂xi


,

(2.3)

where εij is the strain tensor. In elastic soft media, Hooke’s law is satisfied, which provides
the strain-stress relationship as [152]
σij = Cijkl εkl ,

(2.4)

where Cijkl is the compliance matrix which describes the mechanical properties of soft tissue.
The engineering notation for Eq. 2.4 is [152]

 
 σx  C11 C12 C13

 
 σ  C
y

  12 C22 C23

 

 
 σz  C13 C23 C33

=

 
σyz  C14 C24 C34

 

 
σ  C
 zx   15 C25 C35

 
σxy
C16 C26 C36





C14 C15 C16   εx 


 ε 
C24 C25 C26 
y






C34 C35 C36   εz 

,




C44 C45 C46  2εyz 





C45 C55 C56 
 2εzx 


C46 C56 C66
2εxy

(2.5)

where the compliance matrix is a 6 × 6 matrix. For a fully anisotropic medium, the compliance matrix has 21 degrees of freedom. Alternatively, if the medium is isotropic, the
compliance matrix Cijkl contains only 2 degrees of freedom, and the compliance matrix is

15

written as [152]





λ
λ
 σx  λ + 2µ

 
σ   λ
λ + 2µ
λ
 y 

 

 
 σz   λ
λ
λ + 2µ

=

 
σyz   0
0
0

 

 
σ   0
0
0
 xz  

 
σxy
0
0
0





0 0 0   εx 




0 0 0
  εy 




0 0 0   εz 

,


2εyz 
µ 0 0






0 µ 0  2εxz 



0 0 µ
2εxy

(2.6)

where λ and µ are the Lamè constants, also known as the bulk modulus and the shear
modulus, respectively. The tensor notation for Eq. 2.6 is
σij = λεkk δij + 2µεij ,

(2.7)

where δij is the Krönecker symbol, which is written as δij = 1 for i = j and δij = 0 for
i 6= j. Combining Eq. 2.2, 2.3, and 2.7, yields the equation of motion in elastic media [152],
(λ + 2µ)∇(∇ · u) − µ ∇ × (∇ × u) + f = ρ

∂
u,
∂t

(2.8)

where u is the particle displacement, f is the induced body force, and t is the time. Eq. 2.8
is also known as Navier’s equation.
Eq. 2.8 describes elastodynamic waves in a medium without any loss, but the behavior of
soft tissue is not accurately characterized by this model, as tissue is always lossy (viscoelastic)
[79, 116]. To better describe the mechanical response for an acoustic radiation force in
human tissue, a proper loss model is needed. When the Kevin-Voigt model for viscous loss
is included in Eq. 2.8, the Lamè constants λ and µ in Eq. 2.7 are replaced by λ + ηp (∂/∂t)
and µ + ηs (∂/∂t), respectively, where ηp is the compressional viscosity, and ηs is the shear
viscosity [148, 154]. The stress-strain relationship then becomes




∂
∂
σij = λ + ηp
ε δ + 2 µ + ηs
ε .
∂t kk ij
∂t ij

16

(2.9)

After Eq. 2.2, 2.3, and 2.9 are combined, the equation of motion in viscoelastic media is
given by [153]




∂
∂
∂
λ + 2µ + (ηp + 2ηs )
∇(∇ · u) + µ + ηs
∇ × (∇ × u) + f = ρ u.
∂t
∂t
∂t

2.3

(2.10)

Green’s function solutions to the equation of motions

To solve Eqs. 2.8 and 2.10 with Green’s functions that describe the mechanical response to
an applied body force in elastic and viscoelastic media, the displacement u and external body
force f are separated into two parts through Helmholtz decomposition. One part contains
the contributions from scalar potentials that describe the compressional component of the
elastodynamic wave, and the other part contains the contributions from vector potentials that
describe the shear component of the elastodynamic wave [155]. The Helmholtz decomposition
is written as
u = ∇φu + ∇ × ψu ,

(2.11)

f = ∇φf + ∇ × ψf ,

(2.12)

where φu and φf represent the scalar potentials for the displacement field and the body
force, respectively, and ψu and ψf represent the vector potentials for the displacement field
and the body force, respectively.

2.3.1

Green’s function for Navier’s equation in elastic media

For an elastic medium, Navier’s equation in Eq. 2.8 can be separated into two wave equations
that model the scalar and vector components of the wave separately by applying Eqs. 2.11
and 2.12 to Eq. 2.8, which yields
φf
λ + 2µ 2
∂ 2 φu
=
+
∇ φu ,
ρ
ρ
∂t2
ψf
∂ 2 ψu
µ
=
+ ∇2 ψu .
2
ρ
ρ
∂t
17

(2.13)
(2.14)

Eqs. 2.13 and 2.14 describe two distinct wave equations with wave propagation speeds of
p
p
cp = (λ + 2µ)/ρ and cs = µ/ρ, respectively.
By applying a point source to Eq. 2.8, the Green’s function is obtained. A derivation of
this solution is provided in [152]. The Green’s function for Eq. 2.8 is given by
P (r, t) + g S (r, t) + g C (r, t),
gij (r, t) = gij
ij
ij


1
1
r
P
,
gij (r, t) =
γi γj δ0 t −
r
cp
4πρc2p


1
r
1
S
,
gij (r, t) =
(δij − γi γj ) δ0 t −
r
cs
4πρc2s
Z
1
1 r/cs
C
gij (r, t) = −
(3γ γ − δij ) 3
τ δ0 (t − τ ) dτ,
4πρ i j
r r/cp

(2.15)
(2.16)
(2.17)
(2.18)

P (r, t),
where gij (r, t) represents the Green’s function for Navier’s equation in elastic media, gij
S (r, t), and g C (r, t) represent the compressional, shear, and the coupling components of
gij
ij

gij (r, t), respectively, r denotes a spatial point in Cartesian coordinates, i = 1, 2, or 3, which
represents the x, y, and z directions for the displacement, respectively, and j = 1, 2, or
3, which represents the x, y, and z directions for the external body force, respectively. In
the above expressions, r = |r| is the distance between the source point and the observation
point, γi = xi /r, γj = xj /r, δ(·) denotes the Dirac delta function, and H(·) denotes the
Heaviside (unit step) function.
In ultrasound elastography since the ARF push beam is applied along the axial direction,
so the body force is applied in the z direction. The displacement of the propagating shear
wave is also in the z direction. Thus, for the sake of convenience, Eqs. 2.15 - 2.18 are reduced
to the i = 3 and j = 3 case, which represents the shear wave particle displacement in the
z direction induced by the external force applied in the z direction. The simplified Green’s

18

function is then given as
gz (r, t) = gzP (r, t) + gzS (r, t) + gzC (r, t), where


1 z2
r
P
gz (r, t) =
δ t−
,
cp
4πρc2p r3


 
r
1
z2 1
S
δ t−
,
gz (r, t) =
1− 2
r
cs
4πρc2s
r
 2

 



1
3z
1
r
r
C
gz (r, t) =
− 1 3t H t −
−H t−
.
4πρ r2
cp
cs
r

(2.19)
(2.20)
(2.21)
(2.22)

In Eqs. 2.19 - 2.22, gzP (r, t), gzS (r, t), and gzC (r, t) represent the compressional, shear, and the
coupling components, respectively, of the Green’s function for Navier’s equation evaluated
in an elastic medium, where z is the distance between the source point and the observation
point in the z-direction.

2.3.2

Green’s function to the Navier’s equation in viscoelastic media

For a viscoelastic medium, similar operations can be applied to obtain the Green’s function.
Navier’s equation for viscoelastic media in Eq. 2.10 can also be separated into two wave
equations that model the scalar and vector components separately by applying Eqs. 2.11
and 2.12 to Eq. 2.10, which yields


φf
∂ 2 φu
∂
2
=
+ cp + vp
∇2 φu ,
2
ρ
∂t
∂t


2
ψf
∂
∂ ψu
2
=
+ cs + v s
∇2 ψu ,
ρ
∂t
∂t2

(2.23)
(2.24)

where vp = (ηp + 2ηs )/ρ and vs = ηs /ρ.
By applying a point source to Eq. 2.10, an approximate Green’s function is obtained.

19

The approximate Green’s function for Eq. 2.10 is written as
P (r, t) + g S (r, t) + g C (r, t),
gij (r, t) = gij
ij
ij

1
1
1 −
P (r, t) =
p
gij
γi γj e
2
r
4πρcp 2πvp t

(2.25)

(t−r/cp )2 c2p
2vp t
,

(2.26)
2 2

s ) cs
1
1
1 − (t−r/c
2vs t
√
=
(δ
−
γ
γ
)
e
,
ij
i j
r
4πρc2s 2πvs t
C (r, t) = 1 (3γ γ − δ ) 1 I (r, t) + I (r, t) ,
gij
s
ij 3 p
4πρ i j
r


"
p
tc2
(t−r/cp )2 c2p
p
vp t
− 2v
−
t
2vp t
p −e
e
+
Ip (r, t) = √
erf
2
2πcp
S (r, t)
gij

(2.27)
(2.28)
cp t
p
2vp t

!
− erf

cp (t − r/cp )
p
2vp t

!#
,
(2.29)

√
 



(t−r/cs )2 c2s
tc2
s
t
c
t
c
(t
−
r/c
)
vs t  − 2v
−
s
s
s
2vs t
+
s −e
√
Is (r, t) = √
e
erf √
− erf
,
2
2vs t
2vs t
2πcs




(2.30)
where erf(·) is the error function. Eqs. 2.25 - 2.30 are also reduced for the i = 3 and j = 3
case, which represents the particle displacement along the z direction induced by the external

20

force applied along the z-direction. The simplified Green’s function for Eq. 2.10 is given by
gz (r, t) = gzP (r, t) + gzS (r, t) + gzC (r, t), where

(2.31)

(t−r/cp )2 c2p

gzP (r, t) =
gzS (r, t) =
gzC (r, t) =

Ip (r, t) =

1
1
z2 −
2vp t
p
,
e
4πρcp 2πvp t r3


2 2
s ) cs
1
z 2 1 − (t−r/c
1
2vs t
√
1− 2
e
,
4πρcs 2πvs t
r
r
 2


1
3z
1 
−
1
I
(r,
t)
+
I
(r,
t)
,
p
s
4πρ r2
r3





r 2 c2
2
t−
p
p
tcp
cp


νp t
−
−
2νp t
e 2νp − e
 + t erf
√

 2
2πcp

(2.32)
(2.33)
(2.34)
cp t
p
2νp t

!

 

cp t − crp
 ,
− erf  p
2νp t


(2.35)
√
νs t
Is (r, t) = √
2πcs


2
 − tcs
e 2νs





2
t− crs c2s

−
+
2νs t
−e





t
erf
2



c t
√s
2νs t




 
cs t − crs
 ,
− erf  √
2νs t


(2.36)
where gzP (r, t), gzS (r, t), and gzC (r, t) represent the compressional, shear, and the coupling
components, respectively, of the Green’s function for Navier’s equation evaluated in a viscoelastic medium, where z is the distance between the source point and the observation point
in the z-direction.

2.4

Shear wave measurements

In transient ultrasound elastography, shear waves are measured by processing the pulseecho ultrasound signal with a speckle tracking algorithm. In order to measure the local
shear wave particle displacement, multiple frames of RF or IQ data within a 2D FOV are
obtained using pulse-echo ultrasound sequences. For each point in the FOV, the IQ data
within a small spatial window is compared between consecutive frames. The shear wave
particle displacement or velocity is calculated from the RF or IQ data using 1D crosscorrelations [11, 156].
21

2.5

2.5.1

Methods for estimating the mechanical parameters from measured shear waves
The time-of-flight (TOF) method

In shear wave elasticity imaging (SWEI), an estimate of the shear elasticity is often obtained
from the time-of-flight (TOF) of propagating shear waves. In TOF measurements for SWEI,
the arrival times of the shear waves are measured at several successive spatial points located
at the same depth but in different lateral positions. One approach obtains the arrival times
by measuring the time to peak (TTP) [61, 157], which tracks the arrival times of the peaks
for each measured shear wave. Another approach obtains the arrival time delays with crosscorrelations [81, 158, 159]. Of these, the cross-correlation approach is preferred because this
approach is more robust and much less sensitive to noise. The shear wave propagation speed
cs is then calculated from the slope of the linear regression between the lateral locations of
the spatial points and the shear wave arrival times. Another approach that improves the
robustness of the linear regression uses random sample consensus (RANSAC), which updates
the linear fit by randomly generating trial solutions and iteratively updating the model to
eliminate the outliers of the arrival time delays [160, 161]. The shear elasticity µ is then
obtained from µ = ρc2s , where ρ is the density.
Although the TOF method often yields accurate results for SWEI in elastic (lossless)
shear wave phantoms [61, 162–164], this method is based upon some simplifications and
approximations that are responsible for significant bias when applied to shear wave parameter
estimates in soft tissue, which is always lossy (viscoelastic) [79,116]. In viscoelastic media, the
shear wave speed changes during propagation due to shear wave attenuation and dispersion,
so the relationship µ = ρc2s no longer holds [115, 148]. Thus, the shear elasticity estimates
provided by the TOF method are biased by the viscosity. Another problem is that the
TOF method implicitly assumes that the shear wave is a plane wave. However, complicated
diffraction effects are observed in all shear waves induced by a three dimensional (3D) acoustic

22

radiation force [165]. Also, the TOF method yields different values for the estimated shear
elasticity with different transducers, focal depths, and lateral ranges applied to the same
viscoelastic shear wave phantom [166], which suggests that the diffraction of the push beam
and the diffraction of the shear wave also contribute to the bias observed in shear wave
parameter estimates obtained from the TOF method.

2.5.2

The k-space method

Other methods have also been developed for shear wave parameter estimation. One of
these is attenuation measuring ultrasound shearwave elastography (AMUSE) [167]. In the
AMUSE method, a k-space representation of the shear wave is calculated at several successive
spatial points located at the same depth but in different lateral positions using the twodimensional fast Fourier transform (2D-FFT). The k-space representation of the shear wave
data is plotted as function of the angular frequency (ω) and the wavenumber in the xdirection (kx ). Then, the shear wave propagation speed cs is obtained at each frequency by
determining the location of the peak at each frequency in k-space and dividing the value of
ω by the value of kx . The attenuation αs is obtained by measuring the full-width-at-halfmaximum (FWHM) in the kx direction for each angular frequency ω. The AMUSE method
provides a model-free estimate of shear wave speed and attenuation. However, this method
is only effective in the low-frequency range, and diffraction effects are not considered.

23

CHAPTER 3
GPU-BASED GREEN’S FUNCTION SIMULATIONS OF SHEAR WAVES
GENERATED BY AN APPLIED ACOUSTIC RADIATION FORCE IN
ELASTIC AND VISCOELASTIC MEDIA

3.1

Introduction

Shear wave elasticity imaging (SWEI) [42] is an emerging modality for noninvasively
imaging the mechanical properties of soft tissue including breast [81], liver [61], kidney [121],
thyroid [67], heart [168], and blood vessels [150]. In SWEI, an acoustic radiation force push
beam is applied to soft tissue, and the transient shear wave particle displacement is then
measured using high frame-rate ultrasound [42]. Other methods for ultrasound elastography
include acoustic radiation force imaging (ARFI) [41], supersonic shear imaging (SSI) [74],
spatially modulated ultrasound radiation force (SMURF) [86], crawling wave spectroscopy
(CWS) [169, 170], and comb-push ultrasound shear elastography (CUSE) [88].
Numerical simulation approaches are useful for validating and improving shear wave
imaging techniques. At present, most of the numerical simulations of shear waves generated
by an applied acoustic radiation force employ the finite element (FEM) method. FEM
models have been developed for emulating ARFI imaging [138], simulating the shear wave
response to acoustic radiation force in breast [137] and liver [61], modeling the shear wave
propagation induced by an acoustic radiation force in a viscoelastic medium [140], studying
the thermal effects induced by an acoustic radiation force in soft tissue [171], and other
applications [139, 141, 142]. Other simulation options include the finite difference method
[143, 144], the pseudospectral method [146], and the Green’s function approach [147, 148].
The Green’s function approach, which is applicable to shear wave modeling in elastic [147]
and viscoelastic media [148], is appropriate when calculations at only a few observation
points in space are needed, as opposed to other numerical methods that calculate an entire
transient response.
24

Despite the advantages for a small number of points, the Green’s function approach
is very time-consuming for shear wave simulations on a desktop computer. To accelerate
these calculations, a more efficient calculation approach is needed. Since Green’s functions
are readily evaluated in parallel, these calculations are ideal for graphics processing units
(GPUs). The GPU has previously been established as an effective approach for accelerating
medical ultrasound calculations using the angular spectrum method [172], the k-space simulation method [173], and other calculations [172–177]. Similar performance enhancements
are anticipated for shear wave calculations on the GPU.
After the theory of the Green’s function method for the simulation of shear waves generated by an applied acoustic radiation force is introduced, a GPU-based parallel algorithm
that calculates the shear wave motion induced by an acoustic radiation force is described for
simulations in elastic and viscoelastic media. Shear displacements are calculated at several
discrete observation points and also in 2D planes to reveal the effects of shear wave diffraction, attenuation, and dispersion. The computation times of serial and parallel Green’s
function calculations are also evaluated and compared on a desktop computer and on a
high-performance GPU. The convergence of the Green’s function calculations is evaluated
for push beams with different spatial samplings. In addition, simulated shear waves are compared with shear waves measured in different elastic and viscoelastic shear wave phantoms.
The results suggest that the Green’s function approach is ideal for shear wave simulations
with high-performance GPUs.

3.2
3.2.1

Theory
Navier’s Equation

Shear and compressional wave particle displacements induced by an external body force in an
elastic, isotropic, homogeneous, and nearly incompressible model for soft tissue are described

25

by Navier’s equation [152]
(λ + 2µ)∇(∇ · u) − µ ∇ × (∇ × u) + f = ρ

∂
u,
∂t

(3.1)

where u is the local particle displacement, f is the induced body force, and t is the time.
In Eq. 3.1, the bulk modulus and shear modulus are indicated by λ and µ, respectively,
and ρ is the density. When the Kevin-Voigt model for viscous loss is included, Eq. 3.1
becomes [148, 154]




∂
∂
∂
λ + 2µ + (ηp + 2ηs )
∇(∇ · u) + µ + ηs
∇ × (∇ × u) + f = ρ u,
∂t
∂t
∂t

(3.2)

where the parameters ηp and ηs denote the bulk viscosity and shear viscosity, respectively.

3.2.2

Green’s functions calculations

The Green’s function for Navier’s equation evaluated in an elastic medium calculates the
displacement induced by a spatially and temporally impulsive body force. When the acoustic
radiation force push is applied in the z-direction, the z-component of the Green’s function
that describes the resulting displacement in an elastic medium is represented by [152]
gz (r, t) = gzP (r, t) + gzS (r, t) + gzC (r, t), where


r
1 z2
P
δ t−
,
gz (r, t) =
cp
4πρc2p r3

 

1
z2 1
r
S
gz (r, t) =
1− 2
δ t−
,
r
cs
4πρc2s
r

 
 2



1
1
r
3z
r
C
gz (r, t) =
− 1 3t H t −
−H t−
.
4πρ r2
cp
cs
r

(3.3)
(3.4)
(3.5)
(3.6)

In Eqs. 3.3 - 3.6, gzP (r, t), gzS (r, t), and gzC (r, t) represent the compressional, shear, and the
coupling components of the Green’s function for Navier’s equation evaluated in an elastic
medium, respectively, where r denotes a spatial point in Cartesian coordinates, the comp
p
pressional wave speed is cp = (λ + 2µ)/ρ, and the shear wave speed is cs = µ/ρ. In the
above expressions, δ(·) denotes the Dirac delta function, H(·) denotes the Heaviside (unit
26

step) function, r = |r| is the distance between the source point and the observation point,
and z is the distance between the source point and the observation point in the z-direction.
An approximate Green’s function for Navier’s equation evaluated in a viscoelastic medium
is [148]
gz (r, t) = gzP (r, t) + gzS (r, t) + gzC (r, t), where

(3.7)

(t−r/cp )2 c2p

gzP (r, t)

=

gzS (r, t) =
gzC (r, t) =

Ip (r, t) =

1
z2 −
1
2vp t
p
,
e
3
4πρcp 2πvp t r


2 2
s ) cs
1
z 2 1 − (t−r/c
1
2v
t
s
√
1− 2
,
e
4πρcs 2πvs t
r
r

 2

1 
1
3z
−
1
I
(r,
t)
+
I
(r,
t)
,
p
s
4πρ r2
r3





r 2 c2
t−
p
p
tc2
c
p
p
 t
νp t  − 2ν
−
2νp t
p −e
 + erf
e
√
 2
2πcp 

(3.8)
(3.9)
(3.10)
cp t
p
2νp t

!

 

cp t − crp
 ,
− erf  p
2νp t


(3.11)
√
νs t
Is (r, t) = √
2πcs


2
 − tcs
e 2νs





2
t− crs c2s

−
+
2νs t
−e





t
erf
2



c t
√s
2νs t




 
cs t − crs
 ,
− erf  √
2νs t


(3.12)
where νp = (ηp + 2ηs )/ρ, νs = ηs /ρ, and erf(·) denotes the error function.

3.2.3

Simulating shear waves with Green’s functions

To simulate shear waves generated by an acoustic radiation force, the Green’s function and
the acoustic radiation force are convolved in time and space by numerically evaluating the
superposition integral
ZZZ Z
uz (r, t) =
V0

gz (r − r0 , t − τ )f (r0 , τ )dτ dV0 ,

(3.13)

τ

where uz (r, t) is the displacement at r, gz (r, t) is the Green’s function for either the elastic
or viscoelastic model, and f (r0 , t) is the acoustic radiation force. For these calculations, the
27

acoustic radiation force is represented by f (r0 , t) = f (r0 )w(t), where f (r0 ) is the steadystate acoustic radiation force obtained from the intensity calculated in FOCUS, and w(t)
describes the time dependence of the acoustic radiation force. In this model, w(t) is the
same at every point in space, and w(t) is given by w(t) = rect[(t − T /2)/T ]. The expression
rect[(t − T /2)/T ] indicates that the ‘push’ beam is applied starting at time t = 0 and ending
at time t = T , where the variable T describes the duration of the push, and equal weighting
is applied at each point in time.
For shear wave simulations with the elastic model, the convolution in time between the
Green’s function gz (r − r0 , t) in Eqs. 3.3 - 3.6 and w(t) is calculated analytically to eliminate
the delta functions from these calculations. The results are then given by
Z
S
C
hz,T (r, t) = gz (r, t − τ )w(τ )dτ = hP
z,T (r, t) + hz,T (r, t) + hz,T (r, t), where

(3.14)

τ



t − r/cp − T /2
1 z2
rect
,
=
T
4πρc2p r3




t − r/cs − T /2
z2 1
1
S
rect
1− 2
,
hz,T (r, t) =
r
T
4πρc2s
r

 2
1
1
3z
C
−
1
hz,T (r, t) = −
4πρ r2
r3
"
!



 
t − r/cp − T /2
r
1 2 r2
1
2
t − 2 rect
2tT − T H t −
−T
×
+
2
T
2
cp
cp





 
1 r2
r
1
t − r/cs − T /2
2
2
+
− (t − T ) rect
−
2tT − T H t −
,
2 c2s
T
2
cs
hP
z,T (r, t)

(3.15)
(3.16)

(3.17)

where hz,T (r, t) represents the transient displacement of a point body force in the zdirection for the window function w(t) = rect[(t − T /2)/T ] in an elastic medium, and
S
C
hP
z,T (r, t), hz,T (r, t), and hz,T (r, t) are the compressional, shear, and coupling compo-

nents of the resulting transient displacement, respectively. The remaining 3D integral
RRR
hz,T (r − r0 , t)f (r0 )dV0 then describes the spatial convolution between hz,T (r, t) and
V0

the steady-state acoustic radiation force f (r0 ), which is calculated numerically by evaluating
uE
z

Ny Nz
Nx X
X
X

(r, t) =
hz,T (r − r0 , t) f (r0 )∆x∆y∆z.
m=1 n=1 k=1

28

(3.18)

In this triple sum, Nx , Ny , and Nz describe the number of spatial points in each direction
for the simulated steady-state 3D acoustic radiation force, and ∆x, ∆y, and ∆z are the
spatial sampling intervals in the x, y, and z directions for the simulated acoustic radiation
force. The superscript in uE
z (r, t) indicates that this quantity is the shear wave particle
displacement calculated for an elastic medium.
For shear wave simulations with the viscoelastic model, no delta functions appear in
the Green’s function given in Eqs. 3.7 - 3.12, so no analytical simplifications are needed.
For shear wave calculations with Eqs. 3.7 - 3.12 and Eq. 3.13, the convolution in space is
evaluated first, then the convolution in time is computed. The discretized representation of
Eq. 3.13 for numerical calculations in viscoelastic media is


Ny Nz
Nt
Nx X
X
X
X



uVz (r, t) =
gz r − r0 , t − τ f (r0 )∆x∆y∆z  w(τ )∆τ,
i=1

(3.19)

m=1 n=1 k=1

where Nt is the number of time samples in w(t), and ∆τ is the temporal sampling interval.
The superscript in uVz (r, t) indicates that this quantity is the shear wave particle displacement calculated for a viscoelastic medium.

3.3
3.3.1

Numerical Calculations
Acoustic radiation force

In shear wave elasticity imaging (SWEI), shear waves are generated by focused or unfocused
ultrasound beams. For three-dimensional (3D) shear wave calculations in an isotropic, homogeneous simulation model, the pressure field is calculated first, and then the radiation
force is obtained from the computed pressure field, where the radiation force then provides
the input to Eq. 3.1 or Eq. 3.2 for calculations of the shear wave particle displacements.
In each of the following simulations, the pressure, radiation force, and shear wave particle
displacements are all evaluated in 3D.
The 3D pressure fields are calculated in FOCUS, the “Fast Object-oriented C++ Ultrasound Simulator” (http://www.egr.msu.edu/∼fultras-web). In these simulations, the fast
29

nearfield method computes the input pressure plane [130, 131], and then the angular spectrum approach determines the pressure field throughout the remaining 3D volume [132,133].
For each of the shear wave simulations evaluated here, a Philips L7-4 linear ultrasound transducer array (Philips Healthcare, Andover, MA) is modeled, where the array consists of 128
curved transducer elements with an elevation focus of 25 mm. The width of each ultrasound
transducer element is 0.283 mm, the height is 7 mm, and the kerf between adjacent elements
is 0.025 mm. The transducer array face is centered at z = 0 in the x-y plane, the center of
each element is located on the y-axis, and the normal at the center of the array is coincident
with the z-axis.
After the 3D pressure is calculated, the 3D intensity is computed with I = |p|2 /(2ρcp ),
and the acoustic radiation force f (r, t) is determined from [41]
f (r, t) =

2αI(r, t)
,
cp

(3.20)

where α is the ultrasound attenuation coefficient, and I(r, t) is the acoustic intensity.

3.3.2

Parallel implementation

The convolution integral in Eq. 3.13 is very time-consuming when implemented as a serial
calculation. For Nx × Ny × Nz observation points and for Nt time samples, Eq. 3.13 is
evaluated at Nx × Ny × Nz × Nt points in space and time. In practice, this is a problem
for calculations on desktop or laptop computers for large values of Nx × Ny × Nz × Nt .
However, Green’s function calculations are independent at each point in time and at each
point in space, which suggests that these Green’s function calculations are appropriate for
parallel computations. Thus, a graphics processing unit (GPU) is ideal for this problem.
On the GPU, which contains a large number of independent processors, the jobs are broken into several pieces that are distributed to different processors for parallel computations.
When these calculations are implemented properly, GPUs achieve a considerable speedup
compared to a single serial CPU. First, the acoustic radiation force, which is pre-calculated

30

and stored in a data file, is loaded into the CPU (host) memory. Then, the computed acoustic
radiation force is copied into the GPU (device) global memory. The kernel function, which
describes the superposition algorithm, is then called and the job is separated into multiple
threads that are evaluated on different processors in parallel on the GPU. For each thread,
the shear wave particle displacement is calculated at a single observation point and a single moment in time. Each thread obtains the acoustic radiation force from global memory,
whereas the calculated shear wave particle displacement result is first stored temporarily in
local memory for each thread and then copied into global memory after the calculation is
completed. After all of the threads are finished, the results are copied and returned to the
host memory. The CPU performs any post-processing that is needed and then stores the
result in an output data file. These calculations are implemented using the Compute Unified
Device Architecture (CUDA) programming language on an nVidia (Santa Clara, CA) Kepler
K20 GPU.

3.4

Results

Numerical simulations are performed serially on an Intel(R) Core(TM) i7-2600 3.40GHz
CPU and in parallel on an nVidia Tesla K20 GPU. In these simulations, 3D pressure fields
are calculated using the fast nearfield method combined with the angular spectrum approach,
and shear waves are then computed for the 3D acoustic radiation force that is obtained from
the 3D pressure distribution. In each calculation of the acoustic radiation force, the compressional wave speed cp is 1500 m/s, and the attenuation coefficient αp is 0.5 dB/cm/MHz
for the compressional wave component. The excitation frequency for the transducer array
is f0 = 4.09 MHz. For the shear wave simulations, the shear elasticity µs is 2.25 kPa, while
the viscosity ηs is equal to 0 Pa·s for simulations with the elastic model and 1 Pa·s for simulations with the viscoelastic model. The duration of the acoustic radiation force push in
each simulation is 200 µs.

31

0

z (mm)

10
20
30
40
-10

0

10

20

30

40

50

x (mm)
Figure 3.1: The locations of the observation points relative to a push beam focused at 25
mm in the x-z plane with y = 0.
3.4.1

Simulated shear wave particle displacements at discrete observation
points

To demonstrate the similarities and differences between shear waves simulated with the
elastic and viscoelastic models, shear wave particle displacements are calculated at 8 different
observation points located at the same depth. The depth selected for each observation point
is the same as the depth of the focal peak for the push beam. The observation points are
equally spaced in the lateral direction at this depth, where the distance from the peak of the
push beam to each observation point in the x-direction ranges from 5 to 40 mm with a 5 mm
increment in the y = 0 plane. These simulations are evaluated for an f/2 push beam with
40 transducer elements that is focused on-axis at a depth of 25 mm. The focal peak for this
push beam is located at 22.9 mm. The push beam and the locations of the observation points
are shown in Figure 3.1. At each observation point, the shear wave particle displacements
are calculated for 40 ms with a temporal sampling frequency of 10 kHz.
Figure 3.2 plots the simulated shear wave particle displacements for the push beam
shown in Figure 3.1, where the simulated shear wave particle displacements for the elastic

32

1.2
Elastic Medium 3D Push Result
Viscoelastic Medium 3D Push Result

x=5mm

Normalized Displacement

1

0.8

x=10mm
x=15mm
x=20mm
x=25mm
x=30mm
x=35mm
x=40mm

0.6

0.4

0.2

0

−0.2
0

5

10

15

20
t(ms)

25

30

35

40

Figure 3.2: Simulated shear wave particle displacements computed with elastic and viscoelastic models for an applied 3D acoustic radiation force push evaluated at the 8 observation
points indicated in Figure 3.1. The solid lines describe the normalized shear wave particle
displacements calculated for the elastic model with µ = 2.25 kPa and ηs = 0 Pa·s, and the
dashed lines describe the normalized shear wave particle displacements calculated for the
viscoelastic model with µ = 2.25 kPa and ηs = 1 Pa·s.
and viscoelastic models with µ = 2.25 kPa and ηs = 0 or 1 Pa·s are evaluated at the
observation points shown in Figure 3.1. In Figure 3.2, the peaks for different observation
points appear at different times. Figure 3.2 also shows that the peak times for the normalized
shear wave particle displacements calculated with the elastic model are different from the
peak times for the viscoelastic model. Furthermore, at the same observation point, the peaks
calculated for the viscoelastic model arrive later than the peaks for the elastic model. In

33

addition, the differences in the peak times between successive observation points are not
the same for simulations with the elastic model and the viscoelastic model. For shear wave
simulations with the elastic model, the differences in the peak times for adjacent observation
points are 3.34 ms, 3.37 ms, 3.36 ms, 3.35 ms, 3.37 ms, 3.35 ms, and 3.38 ms for the
observation points with x coordinates equal to (5mm and 10mm), (10mm and 15mm), (15mm
and 20mm), (20mm and 25mm), (25mm and 30mm), (30mm and 35mm), and (35mm and
40mm), respectively. For simulations of shear waves with the viscoelastic model, the peak
time differences for the same pairs of observation points are 3.69 ms, 3.61 ms, 3.54 ms,
3.42 ms, 3.33 ms, 3.30 ms, and 3.29 ms, respectively. Thus, the peak time differences are
nearly constant for simulations of shear waves with the elastic model, whereas the peak time
differences for the simulated shear waves with the viscoelastic model decrease monotonically
as the pairs of observation points move further away from the push beam. The amplitudes of
the peak displacements decrease with distance for the elastic model strictly due to the effect
of diffraction, whereas the peak amplitudes for the viscoelastic model decrease due to the
combined effects of diffraction and viscous loss. Also, the peak amplitudes decrease much
more rapidly with the viscoelastic model than with the elastic model.

3.4.2

Computation times for serial and parallel calculations

To evaluate the computation times of the serial and parallel implementations, the simulations are programmed and evaluated in C and CUDA, respectively. Computation times are
determined for shear wave simulations with the elastic and viscoelastic models using the
push beam shown in Figure 3.1. In these simulations, the number of observation points
ranges from 20 to 200 with an increment of 20. The observation points are sampled uniformly between 5 mm and 40 mm in the x direction, the y-coordinate is equal to zero, and
the z coordinate is equal to 22.9 mm. Thus, the end points are the same as in Figure 3.1,
but the region in between is more densely sampled. For each observation point, the shear
wave particle displacements are calculated at 400 temporal samples for 40 ms, which yields

34

a temporal sampling rate of 10 kHz. Eq. 3.13 is evaluated once at each temporal sample
and at each observation point, so the total number of spatial and temporal points at which
the Green’s functions are superposed ranges from 8,000 (20 spatial points × 400 temporal
points) to 80,000 (200 spatial points × 400 temporal points). Furthermore, the discretized
3D acoustic radiation force contains 109 × 55 × 218 = 1,306,910 source points, so the total
number of Green’s function evaluations with Eqs. 3.7 - 3.12 or Eqs. 3.14 - 3.17 ranges from
1.05 × 1010 to 1.05 × 1011 in these simulations. The serial and parallel codes are evaluated
and the computation times are measured and stored. The serial C code is implemented
on a desktop computer with an Intel(R) Core(TM) i7-2600 3.40GHz CPU, and the parallel
CUDA code is implemented on an nVidia Tesla K20 GPU at Michigan State University’s
high performance computing center (https://icer.msu.edu/).
Figure 3.3 describes the computation time versus the number of observation points for
elastic and viscoelastic shear wave simulation models. The results in Figure 3.3 show that the
parallel implementation is much faster than the serial implementation. For the shear wave
simulations with the elastic model, the parallel implementation is about 45 times faster than
the serial implementation, while for shear wave simulations with the viscoelastic model, the
parallel implementation is about 700 times faster than the serial implementation. Moreover,
the simulations of shear waves with the elastic model are much faster than the simulations of
shear waves with the viscoelastic model for both the serial and parallel implementations. In
the serial implementation, the simulation of shear waves with the elastic model is about 50
times faster than the simulation of shear waves with the viscoelastic model. In the parallel
implementation, the simulation of shear waves with the elastic model is around 2 times faster
than the simulation of shear waves with the viscoelastic model.

3.4.3

Convergence of Green’s function calculations

To assess the numerical convergence of these Green’s function calculations, several simulations of shear waves induced by push beams with different spatial sampling rates and

35

106
CUDA code - elastic model
CUDA code - viscoelastic model
C code - elastic model
C code - viscoelastic model

Computation time (s)

105
104
103
102
101
100

0

50

100

150

200

Number of observation points
Figure 3.3: Computation times for serial and parallel implementations of time domain shear
wave calculations with Green’s functions. The solid line with circle markers indicates the
computation times for parallel simulations of shear waves with the elastic model, the dasheddotted line with triangle markers indicates the computation times for parallel simulations
of shear waves with the viscoelastic model, the dashed line with diamond markers indicates
the computation times for serial simulations of shear waves with the elastic model, and the
dotted line with square markers indicates the computation times for serial simulations of
shear waves with the viscoelastic model.
different f-numbers (f/#) are evaluated. Different ultrasound beams are simulated by selectively activating contiguous groups of transducer elements. For unfocused beam calculations,
12 elements are activated with equal phases and amplitudes. For focused beam calculations,
the number of active array elements varies with the f/#. Specifically, 54, 40, and 27 elements
are excited and focused on-axis at a depth of 25 mm for values of the f/# equal to 1.5, 2,
and 3, respectively. Figure 3.4 depicts four different simulated ultrasound beams in the y=0
plane. Figure 3.4a shows the focused ultrasound beam with f/1.5, Figure 3.4b shows the
focused ultrasound beam with f/2, Figure 3.4c shows the focused ultrasound beam with f/3,
and Figure 3.4d shows the unfocused beam (f/∞). The f/2 push beam shown in Figure 3.4b
describes the same configuration depicted in Figure 3.1.
36

0

5

5

5

5

10

10

10

10

15

15

15

15

20

20

z (mm)

0

z (mm)

0

z (mm)

z (mm)

0

20

20

25

25

25

25

30

30

30

30

35

35

35

35

40
-10 -5 0 5
(a)
x (mm)

10

40
-10 -5 0 5
(b)
x (mm)

10

40
-10 -5 0 5
(c)
x (mm)

10

40
-10 -5 0 5
(d)
x (mm)

10

Figure 3.4: 2D cross-sections of 3D ultrasound beams simulated in FOCUS. Figures (3.4a),
(3.4b), (3.4c), and (3.4d) show the push beams with f/1.5, f/2, f/3, and f/∞, respectively.
In (3.4a)-(3.4d), each beam is normalized by the maximum intensity value.
For each push beam, the acoustic radiation force is calculated with spatial samplings of
4λ, 2λ, λ, λ/2, λ/4, λ/8, λ/16, and λ/18, where λ is the wavelength of the compressional
wave (λ = cp /f0 = 0.3667 mm). The push beams with spatial samplings less than or equal
to λ/2 are calculated with FOCUS as described in Section 3.3, and the push beams that are
more coarsely sampled are obtained by down-sampling the push beam that is computed with
λ/2 sampling. The push beam sampled at 4λ, which is the most coarsely sampled push beam,
consists of Nxs × Nys × Nzs = 13 × 7 × 27 = 2,457 source points, while the most finely sampled
push beam sampled at λ/18 consists of Nxs × Nys × Nzs = 981 × 491 × 1963 = 945,520,173
source points. The simulated shear wave particle displacement induced by the λ/18 sampled
push beam is used as the reference, and the simulated shear wave particle displacements
induced by all of the other push beams are compared with the reference shear wave particle
displacement to evaluate the convergence of these Green’s function calculations. The percent
difference between the two shear wave particle displacement waveforms is defined by
% difference =

||u(r, t) − uref (r, t)||
,
||uref (r, t)||

(3.21)

where u(r, t) is the simulated shear wave particle displacement induced by a push beam at a
single observation point, uref (r, t) is the reference shear wave particle displacement induced
37

18
Unfocused push beam
Focused push beam with f/# = 3
Focused push beam with f/# = 2
Focused push beam with f/# = 1.5

16

Percent difference

14
12
10
8
6
4
2
0
/2

/4

/8

/16

Spatial sampling of the push beam
Figure 3.5: Percent differences calculated for shear wave particle displacements simulated
with the elastic model. The shear elasticity and shear viscosity values in these simulations
are µ = 2.25 kPa and ηs = 0 Pa·s, where the push beams are calculated with different
spatial samplings. The solid line with circle markers indicates the percent difference for
the shear wave particle displacements induced by an unfocused push beam with f/∞, the
dotted line with square markers indicates the percent difference for the shear wave particle
displacements induced by a focused push beam with f/3, the dash-dotted line with triangle
markers indicates the percent difference for the shear wave particle displacements induced
by a focused push beam with f/2, and the dashed line with diamond markers indicates the
percent difference for the shear wave particle displacements induced by a focused push beam
with f/1.5.
by the push beam sampled at λ/18 at the same observation point, and || · || denotes the
L2 -norm of the time-domain waveform.
For the shear wave particle displacement simulated by each push beam, the percent
difference is evaluated at 8 observation points. For each push beam, the observation points

38

are all located at the same depth, and the depth selected for each observation point is the
same as the depth of the focal peak. For the push beams with f/1.5, f/2, f/3, and f/∞,
these depths are equal to 24.5 mm, 22.9 mm, 20.0 mm, and 16.6 mm, respectively. The
observation points are equally spaced in the lateral direction, where the distance from the
peak of the push beam to each observation point in the x-direction ranges from 5 to 40 mm
with a 5 mm increment in the y = 0 plane. The largest percent difference is retained to
represent the worst case value across these 8 locations. The percent difference is evaluated
for simulated shear waves with the elastic model using µs = 2.25 kPa and ηs = 0 Pa·s
and for simulated shear waves with the viscoelastic model using µs = 2.25 kPa and ηs = 1
Pa·s for the four push beams with f/1.5, f/2, f/3, and f/∞ described in section 3.3.1. The
results are illustrated in Figure 3.5, which shows the percent difference in the shear wave
particle displacements calculated using the elastic model for the four push beams evaluated
with different spatial samplings. In each of these simulations, the percent difference decreases
when the spatial sampling rate increases. For the f/1.5 push beams sampled at λ, the percent
difference is 16.4%, the percent difference is 7.04% when the spatial sampling is λ/2, and the
percent difference is 2.50% when the spatial sampling of the push beam is λ/16. The percent
differences for the f/2 push beams sampled at λ, λ/2, and λ/16 are 13.5%, 4.53%, and 0.98%,
respectively. The percent differences for the f/3 push beams sampled at λ, λ/2, and λ/16
are 9.96%, 3.63%, and 0.23%, respectively. The percent differences for the unfocused push
beams sampled at 4λ, λ/2, and λ/16 are 6.00%, 2.70%, and 0.29%, respectively.
The percent differences in the shear wave particle displacements simulated for the viscoelastic model with µ = 2.25 kPa and ηs = 1 Pa·s are displayed in Figure 3.6 for the same
four push beams. Similar to the Green’s function simulations of shear waves with the elastic
model, the percent difference decreases when the spatial sampling increases. For the f/1.5
push beam sampled at λ, the percent difference is 6.15%, the percent difference is 4.74%
when the spatial sampling is λ/2, and the percent difference is 1.00% when the spatial sampling of the push beam is λ/16. The percent differences for the f/2 push beam sampled at

39

10
Unfocused push beam
Focused push beam with f/# = 3
Focused push beam with f/# = 2
Focused push beam with f/# = 1.5

Percent difference

8

6

4

2

0
/2

/4

/8

/16

Spatial sampling of the push beam
Figure 3.6: Percent differences calculated for shear wave particle displacements simulated
with the elastic model. The shear elasticity and shear viscosity values in these simulations
are µ = 2.25 kPa and ηs = 1 Pa·s, where the push beams are calculated with different spatial
samplings. The descriptions of the markers and line types are the same as those given in the
caption for Figure 3.5.
λ, λ/2, and λ/16 are 2.56%, 1.26%, and 0.40%, respectively. The percent differences for the
f/3 push beam sampled at λ, λ/2, and λ/16 are 1.15%, 0.92%, and 0.19%, respectively. The
percent differences for the unfocused push beam sampled at λ, λ/2, and λ/16 are 1.35%,
0.83%, and 0.34%, respectively.

40

3.4.4

Comparisons between simulated and measured shear waves

To validate the viscoelastic model, shear waves are measured in viscoelastic shear wave
phantoms and compared with the simulated shear waves. Shear waves were measured with a
Verasonics ultrasound system (Verasonics Inc., Redmond, WA) and an L7-4 linear ultrasound
transducer array, which is the same transducer array that is described and modeled in section
3.3.1. To induce shear waves, an f/2 push beam focused at 25 mm with a center frequency
of 4.09 MHz is generated when 40 contiguous elements of the transducer array are activated
with the same f/2 focused beam shown in Figure 3.1. The push beam is applied for 200
µs, and then the Verasonics system is switched to plane wave imaging mode for shear wave
measurements. In-phase/quadrature (IQ) data is measured for 10 ms with a frame rate of
11.76 kHz in a two-dimensional (2D) region ranging from -19.40 mm to 19.40 mm in the
lateral direction and from 0 to 39.58 mm in the axial direction with a spatial sampling of 0.154
mm in each direction. Then, shear wave particle velocities are obtained from the IQ data by
calculating one-dimensional (1D) autocorrelations between consecutive frames [11]. A 2D
median filter with a 3 by 3 kernel is then applied to the measured shear wave particle velocities
for each frame. Three different viscoelastic shear wave phantoms from Computerized Imaging
Reference Systems (CIRS), Inc (Norfolk, Virginia, US) are used for these measurements,
namely shear wave phantoms E3962-1, E3962-2, and E3962-3.
After the shear waves in the viscoelastic shear wave phantoms are measured, the results
are compared to simulated shear waves evaluated within the same 2D region at three different
points in time. For phantoms E3962-1 and E3962-2, shear wave particle displacements
are displayed at 3 ms, 6 ms, and 9 ms, while for phantom E3962-3, shear wave particle
displacements are displayed at 2 ms, 3 ms, and 4 ms. Different time points are required
for phantom E3962-3 because the shear elasticity for this phantom is much larger, so the
shear waves propagate out of the FOV much more quickly. The shear waves are simulated for
several different combinations of the shear elasticity and shear viscosity. For phantom E39621, the shear elasticity and the shear viscosity are relatively small, so the shear elasticity in the

41

simulations ranges from 0 to 0.5 kPa with an increment of 0.1 kPa, and the shear viscosity
rages from 0 to 0.5 Pa·s with an increment of 0.1 Pa·s. For shear wave phantoms E3962-2
and E3962-3, the shear elasticity in the simulations ranges from 2 to 5 kPa with an increment
of 0.1 kPa, and the shear viscosity rages from 0 to 2.5 Pa·s with an increment of 0.1 Pa·s.
For each combination of parameters, the simulation of the shear wave particle velocity in a
2D plane at one point in time takes 25 s on an nVidia Tesla K20 GPU. After the shear wave
particle velocities are simulated on the GPU, the peak values of the measured and simulated
shear wave particle velocities are both normalized to 1, and the percent difference between
the simulated and measured shear wave particle velocities is evaluated with
% difference =

||usim (r, t) − umeas (r, t)||
,
||umeas (r, t)||

(3.22)

where usim (r, t) represents the simulated shear wave particle velocity in a 2D plane at time t,
umeas (r, t) represents the measured shear wave particle velocity in the same 2D plane at time
t, and || · || denotes the L2 -norm. The percent difference is calculated for each combination
of the shear elasticity and shear viscosity parameters, and the combination of shear wave
parameters that provides the minimum percent difference is determined from the results.
Figure 3.7 shows the measured shear wave particle velocities in phantom E3962-1 and the
corresponding simulated shear wave particle velocities that minimize the percent difference
in Eq. 3.22. For all three points in time, the material parameters that provide the minimum
percent difference are µ = 0.3 kPa and ηs = 0.1 Pa·s. In Figure 3.7(a), the location of the
peak shear wave particle velocity in the right half of the image is (x, z) = (2.15 mm, 22.18
mm), while in Figure 3.7(d), the location of the peak shear wave particle velocity in the right
half of the image is (x, z) = (2.00 mm, 23.56 mm). Also, in Figure 3.7(a), the full width
at half maximum (FWHM) of the shear wave in the x direction evaluated at the peak on
the right is 1.08 mm, and in Figure 3.7(d), the FWHM of the shear wave in the x direction
evaluated at the peak on the right is 1.14 mm. In Figure 3.7(b), the location of the peak
shear wave particle velocity in the right half of the image is (x, z) = (3.69 mm, 22.79 mm),
while in Figure 3.7(e), the maximum location of the peak shear wave particle velocity in the
42

time = 3 ms

5

15
20
25
-10

0
10
x(mm)

20

15
20

30
-20
(b)

15
20
25

-10

0
10
x(mm)

20

30
-20
(c)

5

5

10

10

10

20
25

30
-20
(d)

z(mm)

5

15

15
20
25

-10

0
10
x(mm)

20

30
-20
(e)

time = 9 ms

10

25

z(mm)

z(mm)

30
-20
(a)

5

10
z(mm)

z(mm)

10

time = 6 ms

z(mm)

5

-10

0
10
x(mm)

20

-10

0
10
x(mm)

20

15
20
25

-10

0
10
x(mm)

20

30
-20
(f)

Figure 3.7: Measured shear wave particle velocities in phantom E3962-1 and the corresponding simulated shear wave particle velocities that minimize Eq. 3.22 at three different times.
Figs. 3.7(a), 3.7(b), and 3.7(c) show the shear wave particle velocities measured at 3 ms, 6
ms, and 9 ms, respectively. Figs. 3.7(d), 3.7(e), and 3.7(f) show the corresponding simulated
shear wave particle velocities computed with µ = 0.3 kPa and ηs = 0.1 Pa·s at 3 ms, 6 ms,
and 9 ms, respectively.
right half of the image is (x, z) = (3.85 mm, 22.95 mm). Also, in Figure 3.7(b), the FWHM
of the shear wave in the x direction evaluated at the peak on the right is 1.55 mm, and in
Figure 3.7(e), the FWHM of the shear wave in the x direction evaluated at the peak on the
right is 1.57 mm. In Figure 3.7(c), the location of the peak shear wave particle velocity in the
right half of the image is (x, z) = (5.24 mm, 21.25 mm), while in Figure 3.7(f), the location
of the peak shear wave particle velocity in the right half of the image is (x, z) = (5.54 mm,
22.48 mm). Also, in Figure 3.7(c), the FWHM of the shear wave in the x direction evaluated
at the peak on the right is 1.98 mm, and in Figure 3.7(f), the FWHM of the shear wave
in the x direction evaluated at the peak on the right is 1.91 mm. In each pair of images,

43

time = 3 ms

5

15
20
25
-10

0
10
x(mm)

20

15
20

30
-20
(b)

15
20
25

-10

0
10
x(mm)

20

30
-20
(c)

5

5

10

10

10

20
25

30
-20
(d)

z(mm)

5

15

15
20
25

-10

0
10
x(mm)

20

30
-20
(e)

time = 9 ms

10

25

z(mm)

z(mm)

30
-20
(a)

5

10
z(mm)

z(mm)

10

time = 6 ms

z(mm)

5

-10

0
10
x(mm)

20

-10

0
10
x(mm)

20

15
20
25

-10

0
10
x(mm)

20

30
-20
(f)

Figure 3.8: Measured shear wave particle velocities in phantom E3962-2 and the corresponding simulated shear wave particle velocities that minimize Eq. 3.22 at three different times.
Figs. 3.8(a), 3.8(b), and 3.8(c) show the shear wave particle velocities measured at 3 ms, 6
ms, and 9 ms, respectively. Figs. 3.8(d), 3.8(e), and 3.8(f) show the corresponding simulated
shear wave particle velocities with µ = 2.3 kPa and ηs = 0.1 Pa·s at 3 ms, 6 ms, and 9 ms,
respectively.
the peak locations and the FWHM values at these peak locations are similar. The images
showing the measured (Figure 3.7a - 3.7c) and simulated (Figure 3.7d - 3.7f) shear wave
particle velocities are also visually similar for phantom E3962-1.
Figure 3.8 shows the measured shear wave particle velocities in phantom E3962-2 and the
corresponding simulated shear wave particle velocities that minimize the percent difference
in Eq. 3.22. For all three points in time, the material parameters that provide the minimum
percent difference are µ = 2.3 kPa and ηs = 0.1 Pa·s. In Figure 3.8(a), the location of the
peak shear wave particle velocity in the right half of the image is (x, z) = (4.93 mm, 22.33
mm), while in Figure 3.8(d), the location of the peak shear wave particle velocity in the right

44

half of the image is (x, z) = (5.08 mm, 23.41 mm). Also, in Figure 3.8(a), the full width
at half maximum (FWHM) of the shear wave in the x direction evaluated at the peak on
the right is 1.23 mm, and in Figure 3.8(d), the FWHM of the shear wave in the x direction
evaluated at the peak on the right is 1.17 mm. In Figure 3.8(b), the location of the peak
shear wave particle velocity in the right half of the image is (x, z) = (9.70 mm, 23.56 mm),
while in Figure 3.8(e), the location of the peak shear wave particle velocity in the right half of
the image is (x, z) = (9.55 mm, 23.10 mm). Also, in Figure 3.8(b), the FWHM of the shear
wave in the x direction evaluated at the peak on the right is 1.62 mm, and in Figure 3.8(e),
the FWHM of the shear wave in the x direction evaluated at the peak on the right is 1.59
mm. In Figure 3.8(c), the location of the peak shear wave particle velocity in the right half
of the image is (x, z) = (14.17 mm, 21.71 mm), while in Figure 3.8(f), the location of the
peak shear wave particle velocity in the right half of the image is (x, z) = (14.32 mm, 22.48
mm). Also, in Figure 3.8(c), the FWHM of the shear wave in the x direction evaluated at
the peak on the right is 1.80 mm, and in Figure 3.8(f), the FWHM of the shear wave in the
x direction evaluated at the peak on the right is 1.94 mm. In each pair of images, the peak
locations and the FWHM values at these peak locations are similar. The images showing
the measured (Figure 3.8a - 3.8c) and simulated (Figure 3.8d - 3.8f) shear wave particle
velocities are also visually similar for phantom E3962-2.
Figure 3.9 shows the measured shear wave particle velocities in phantom E3962-3 and the
corresponding simulated shear wave particle velocities that minimize the percent difference
in Eq. 3.22. For the shear waves that are measured and simulated at 2 ms, the material
parameters that provide the minimum percent difference are µ = 4.9 kPa and ηs = 1.5
Pa·s. For the shear waves that are measured and simulated at 3 ms and 4 ms, the material
parameters that provide the minimum percent difference are µ = 4.5 kPa and ηs = 1.5
Pa·s. In Figure 3.9(a), the location of the peak shear wave particle velocity in the right half
of the image is (x, z) = (6.00 mm, 19.25 mm), while in Figure 3.9(d), the location of the
peak shear wave particle velocity in the right half of the image is (x, z) = (5.85 mm, 19.40

45

time = 2 ms

5

15
20
25
-10

0
10
x(mm)

20

15
20

30
-20
(b)

15
20
25

-10

0
10
x(mm)

20

30
-20
(c)

5

5

10

10

10

20
25

30
-20
(d)

z(mm)

5

15

15
20
25

-10

0
10
x(mm)

20

30
-20
(e)

time = 4 ms

10

25

z(mm)

z(mm)

30
-20
(a)

5

10
z(mm)

z(mm)

10

time = 3 ms

z(mm)

5

-10

0
10
x(mm)

20

-10

0
10
x(mm)

20

15
20
25

-10

0
10
x(mm)

20

30
-20
(f)

Figure 3.9: Measured shear wave particle velocities in phantom E3962-3 and the corresponding simulated shear wave particle velocities that minimize Eq. 3.22 evaluated at three different
times. Figs. 3.9(a), 3.9(b), and 3.9(c) show the shear wave particle velocities measured at
2 ms, 4 ms, and 6 ms, respectively. Figs. 3.9(d), 3.9(e), and 3.9(f) show the corresponding
simulated shear wave particle velocities computed with µ = 4.9 kPa and ηs = 1.5 Pa·s at 2
ms, and computed with µ = 4.5 kPa and ηs = 1.5 Pa·s at 3 ms and 4 ms, respectively.
mm). Also, in Figure 3.9(a), the full width at half maximum (FWHM) of the shear wave
in the x direction evaluated at the peak on the right is 3.04 mm, and in Figure 3.9(d), the
FWHM of the shear wave in the x direction evaluated at the peak on the right is 3.37 mm.
In Figure 3.9(b), the location of the peak shear wave particle velocity in the right half of
the image is (x, z) = (7.39 mm, 17.25 mm), while in Figure 3.9(e), the location of the peak
shear wave particle velocity in the right half of the image is (x, z) = (8.01 mm, 17.71 mm).
Also, in Figure 3.9(b), the FWHM of the shear wave in the x direction evaluated at the
peak on the right is 3.55 mm, and in Figure 3.9(e), the FWHM of the shear wave in the
x direction evaluated at the peak on the right is 4.24 mm. In Figure 3.9(c), the location

46

of the peak shear wave particle velocity in the right half of the image is (x, z) = (10.93
mm, 15.25 mm), while in Figure 3.9(f), the location of the peak shear wave particle velocity
in the right half of the image is (x, z) = (10.47 mm, 15.86 mm). Also, in Figure 3.9(c),
the FWHM of the shear wave in the x direction evaluated at the peak on the right is 4.62
mm, and in Figure 3.9(f), the FWHM of the shear wave in the x direction evaluated at the
peak on the right is 4.81 mm. In each pair of images, the peak locations are quite similar,
but the FWHM values at these peak locations differ somewhat, since the FWHM for the
measured shear wave is consistently smaller than the FWHM for the simulated shear wave.
The above results indicate that, although some differences are observed in the FWHM values
for the shear wave particle velocities, the peak locations are similar, and these differences and
similarities are also observed in visual comparisons between Figures 3.9(a-c) and 3.9(d-f).

3.5
3.5.1

Discussion
Diffraction, attenuation, and dispersion

The results shown in Figure 3.2 reveal the complicated interaction between diffraction, attenuation, and dispersion in the shear wave particle displacements induced by a 3D acoustic
radiation force. Figure 3.2 shows that the peak amplitude of the shear wave particle displacement decreases as the observation point moves away from the push beam. For simulations
of shear waves with the elastic model, the amplitude of the shear wave particle displacement
decreases solely due to diffraction. The peak time differences between successive observation
points are approximately the same, which indicates that the propagation speed of the shear
wave remains constant along the line of observation points shown in Figure 3.1. Moreover,
the values of the full width at half maximum (FWHM) for the shear wave particle displacement calculated at the observation points shown in Figure 3.1 are 1.06 ms, 1.02 ms, 0.99
ms, 0.97 ms, 1.00 ms, 1.02 ms, 1.06 ms, and 1.08 ms for the pairs of observation points
with x coordinates equal to 5 mm, 10 mm, 15 mm, 20 mm, 25 mm, 30 mm, 35 mm, and
40 mm, respectively. Thus, the values of the FWHM are similar at each of these observa-

47

tion points, which suggests that there is no significant waveform spreading due to dispersion
along this line of points for simulated shear wave particle displacements evaluated with the
elastic model. For simulations of shear waves with the viscoelastic model, the amplitude of
the shear wave particle displacement decreases much faster than with the elastic model, as
the amplitude of the shear wave particle displacement decreases due to the combined effects
of diffraction and attenuation. Furthermore, the differences between the times at which
the peaks occur for successive observation points are monotonically decreasing as the observation points move away from the push beam. This suggests that the propagation speed
of the shear wave particle displacement decreases as the shear wave propagates due to the
frequency dependence of the phase velocity. Moreover, the values of the full width at half
maximum (FWHM) for the shear wave particle displacements calculated at the observation
points shown in Figure 3.1 are 5.50 ms, 6.70 ms, 7.33 ms, 7.68 ms, 7.94 ms, 8.22 ms, 8.54 ms,
and 8.89 ms for the observation points with x coordinates equal to 5 mm, 10 mm, 15 mm, 20
mm, 25 mm, 30 mm, 35 mm, and 40 mm, respectively. These values for the FWHM, which
increase monotonically, also demonstrate that the shear wave particle velocity spreads out
in time as the shear wave propagates, thereby providing additional evidence of shear wave
dispersion in the viscoelastic model.

3.5.2

Serial vs parallel implementations

The computation times for the serial and parallel calculations shown in Figure 3.3 indicate
that the GPU accelerates these Green’s function calculations by factors of 45 and 700 for
simulations of shear waves with the elastic and viscoelastic models, respectively. Thus, the
GPU achieves a significant reduction in the computation time. For instance, the computation
time for the shear waves in a 252 by 257 plane shown in Figs. 3.7-3.9 is around 26.3 seconds
on an nVidia Tesla K20 GPU, whereas the serial implementation of the same calculation
takes 5.27 hours on a desktop computer with an Intel(R) Core(TM) i7-2600 3.40GHz CPU.
This suggests that these Green’s function calculations are preferably performed on a high-

48

performance GPU. The results shown in Figure 3.3 also reveal that both serial and parallel
calculations are faster for simulations of shear waves with the elastic model because the
calculation of the error function in Eqs. 3.7 - 3.12 is much more time consuming than the
calculation of the Heaviside function in Eqs. 3.14 - 3.17.

3.5.3

Convergence

Figures 3.5 and 3.6 show that the percent difference is larger for shear wave particle displacements calculated using push beams with lower f/#. Smaller percent differences are achieved
in simulations of arrays with larger f/# because for push beams with lower f/#, the energy of the acoustic radiation force is more concentrated in the focal region, so finer spatial
samplings are required to capture the high frequency components in the spatial domain.
Figures 3.5 and 3.6 also indicate that when the spatial sampling of the push beam is
greater than or equal to λ/2, the percent difference for the simulated shear wave particle
displacements with the elastic model is usually somewhat higher than for simulations of
shear waves with the viscoelastic model. However, for push beams sampled more finely than
λ/2, the percent differences are similar for shear wave simulations evaluated with the elastic
and viscoelastic models.
For most of the results shown in Figures 3.5 and 3.6, the percent difference decreases
monotonically as the spatial sampling decreases. Two exceptions to this trend are observed
in Figure 3.6, where the percent difference for push beams sampled at λ/2 is smaller than the
percent difference for push beams sampled at λ/4 and λ/8. Similar behavior is also observed
in the convergence of the fast nearfield method for rectangular and circular pistons [130,131],
where the overall trend is that the difference relative to the references is decreasing, but local
minima are also observed in several locations. To determine whether the local minima in
Figure 3.6 are unique to the particular parameter combinations evaluated here, additional
simulations of shear waves in viscoelastic models were performed using f/1.5 push beams
evaluated with different values of the shear elasticity and shear viscosity. Four different

49

parameter combinations were evaluated: µs = 2.25 kPa and ηs = 2 Pa·s, µs = 2.25 kPa and
ηs = 0.1 Pa·s, µs = 6 kPa and ηs = 1 Pa·s, and µs = 1 kPa and ηs = 1 Pa·s. For each
parameter combination, the percent difference for push beams with spatial samplings of mλ
with m = 1 to 4 with an increment of 1 and λ/n with n = 2 to 18 with an increment of 1
were evaluated. In all of these simulations, a local minimum for the percent difference was
also observed for the push beam sampled at λ/2. In particular, the percent difference for
the shear wave particle displacement simulated using the push beam sampled at λ/2 was
consistently smaller than the push beam sampled from 4λ to λ and from λ/3 to λ/8 in each
of these results. This suggests that λ/2 is a reasonable default value for the spatial sampling,
especially for push beams with smaller f/# evaluated with the viscoelastic model. For shear
wave simulations with the elastic model, the percent difference is below 5% for the f/2, f/3,
and f/∞ push beams when the spatial sampling is λ/2, and the percent difference is below
5% for the f/1.5 push beam when the spatial sampling is λ/4. For shear wave simulations
evaluated with the viscoelastic model, the percent difference is below 5% for the f/2, f/3, and
f/∞ push beams when the spatial sampling is 2λ, and the percent difference is below 5% for
the push beam with f/1.5 when the spatial sampling is λ/2. This information is important
for simulations of shear waves because the accuracy of the shear wave particle displacement
depends on the sampling of the push beam.
If a smaller percent difference is desired, e.g., 1%, then a much finer spatial sampling is
required. For shear wave simulations evaluated with the elastic model, the percent difference
is greater than 1% for the f/1.5 push beam when the spatial sampling is greater than or equal
to λ/16, the percent difference is below 1% for the f/2 push beam when the spatial sampling
is λ/16, the percent difference is below 1% for the f/3 push beam when the spatial sampling
is λ/8, and the percent difference is below 1% for the f/∞ push beam when the spatial
sampling is equal to λ/4. For shear wave simulations evaluated with the viscoelastic model,
the percent difference is 1% for the f/1.5 push beam when the spatial sampling is λ/16, the
percent difference is below 1% for the f/2 push beam when the spatial sampling is λ/16,

50

FWHM (mm) at 3ms
FWHM (mm) at 6ms
FWHM (mm) at 9ms

Phantom E3962-1
Measured Simulated
1.08
1.14
1.55
1.57
1.98
1.91

Phantom E3962-2
Measured Simulated
1.23
1.17
1.62
1.59
1.80
1.94

Table 3.1: The FWHM values of the measured and simulated shear wave particle velocities
in the x direction evaluated at the peak at 3 ms, 6 ms, and 9 ms for phantom E3962-1 and
E3962-2.
and the percent difference is below 1% for the f/3 and f/∞ push beams when the spatial
sampling is λ/2.

3.5.4

Comparisons between measured and simulated shear waves

Fig. 3.7 shows the measured and simulated shear waves for phantom E3962-1. The lateral
positions (x) of the peak shear wave particle velocities at 3 ms, 6 ms, and 9 ms are 2.15 mm,
3.85 mm, and 5.24 mm, respectively, for the measured shear waves, and are 2.00 mm, 3.85
mm, and 5.54 mm, respectively, for the simulated shear waves. The depths (z) of the peak
shear wave particle velocities at 3 ms, 6 ms, and 9 ms are 22.18 mm, 22.79 mm, and 21.25
mm, respectively, for the measured shear waves, and are 23.56 mm, 22.95 mm, and 22.48
mm, respectively, for the simulated shear waves. These values indicate that the peak shear
wave particle velocity remains at about the same depth during propagation. In addition,
as shown in Table 3.1, The FWHM values of the shear wave in the x direction evaluated
at the peak at 3 ms, 6 ms, and 9 ms are 1.08 mm, 1.55 mm, and 1.98 mm, respectively,
for the measured shear waves, and are 1.14 mm, 1.57 mm, and 1.91 mm, respectively, for
the simulated shear waves. These results demonstrate that the FWHM values increase by a
significant amount as the shear wave propagates due to shear wave diffraction and waveform
spreading despite the relatively small value for the shear viscosity.
Fig. 3.8 shows the measured and simulated shear waves for phantom E3962-2. The lateral
positions (x) of the peak shear wave particle velocities at 3 ms, 6 ms, and 9 ms are 4.93
mm, 9.70 mm, and 14.17 mm, respectively, for the measured shear waves, and are 5.08 mm,
51

FWHM (mm) at 2ms
FWHM (mm) at 3ms
FWHM (mm) at 4ms

Phantom E3962-3
Measured Simulated
3.04
3.37
3.55
4.24
4.62
4.81

Table 3.2: The FWHM values of the measured and simulated shear wave particle velocities
in the x direction evaluated at the peak at 3 ms, 6 ms, and 9 ms for phantom E3962-3.
9.55 mm, and 14.32 mm, respectively, for the simulated shear waves. The depths (z) of
the peak shear wave particle velocities at 3 ms, 6 ms, and 9 ms are 22.33 mm, 23.56 mm,
and 21.71 mm, respectively, for the measured shear waves, and are 23.40 mm, 23.10 mm,
and 22.48 mm, respectively, for the simulated shear waves. These results indicate that the
peak shear wave particle velocity remains at about the same depth during propagation. The
above results also indicate that the positions of the peak shear wave particle velocity for the
measured and simulated shear waves are similar. In addition, as shown in Table 3.1, the
FWHM values for the shear wave particle velocities in the x direction for the peak on the
right evaluated at 3 ms, 6 ms, and 9 ms are 1.23 mm, 1.62 mm, and 1.80 mm, respectively,
for the measured shear waves, and are 1.17 mm, 1.59 mm, and 1.94 mm, respectively, for
the simulated shear waves. These results demonstrate that the FWHM values increase by a
significant amount as the shear wave propagates due to shear wave diffraction and waveform
spreading despite the relatively small value for the shear viscosity.
Fig. 3.9 shows the measured and simulated shear waves for phantom E3962-3. The lateral
positions (x) of the peak shear wave particle velocities at 2 ms, 3 ms, and 4 ms are 6.00 mm,
7.39 mm, and 10.93 mm, respectively, for the measured shear waves, and are 5.85 mm, 8.01
mm, and 10.47 mm, respectively, for the simulated shear waves. The depths (z) of the peak
shear wave particle velocities at 2 ms, 3 ms, and 4 ms are 19.25 mm, 17.25 mm, and 15.25
mm, respectively, for the measured shear waves, and are 19.40 mm, 17.71 mm, and 15.86
mm, respectively, for the simulated shear waves. These results indicate that the peak shear
wave particle velocity is moving upward in Fig. 3.9 during propagation. The above results
also indicate that the locations for the peak shear wave particle velocities of measured and
52

simulated shear waves have similar positions. In addition, as shown in Table 3.2, the FWHM
of the shear wave in the x direction evaluated for the peak value at 2 ms, 3 ms, and 4 ms are
3.04 mm, 3.55 mm, and 4.62 mm, respectively, for the measured shear waves, and are 3.37
mm, 4.24 mm, and 4.81 mm, respectively, for the simulated shear waves. These results show
that, for larger values of the shear viscosity, there is a significant difference in the FWHM
values. Possible explanations for this are that a different model is required (e.g. a fractional
Voigt model) for larger shear viscosity values and/or that a different objective function that
emphasizes agreement in the FWHM at the expense of diminished agreement at the peak
location is needed. The resolution of these and other related issues is a topic of ongoing
research.

3.5.5

Intended applications

These GPU-accelerated shear wave simulations are intended for several different types of calculations. For example, these shear wave simulations are useful for evaluating different push
beams and different material parameter combinations. These simulations are also applicable
to parametric characterizations of various methods for estimating shear wave parameters,
including the time-of-flight approach [61,81,157–159] and the k-space approach [167]. In addition, these routines can be included in new model-based approachs to minimize the difference
between measured and simulated shear waves to obtain more accurate estimates of shear wave
parameters. Hardware acceleration is essential for the optimization and for repeated shear
wave calculations, because these calculations are otherwise very time-consuming. These calculations could potentially be incorporated into FOCUS, and we also anticipate that other
applications will be developed once these routines are made available.

3.6

Conclusion

Shear waves induced by an acoustic radiation force are evaluated using Green’s functions
for simulations with elastic and viscoelastic models. The 3D acoustic radiation force is first

53

simulated using the fast nearfield method and the angular spectrum approach in FOCUS,
and then shear waves are calculated by superposing the contributions from all of the source
points in the 3D push beam using Green’s functions. Next, these simulations of shear waves
in elastic and viscoelastic models are accelerated using a high-performance GPU. Diffraction
effects are observed in the simulated shear wave particle displacements in both elastic and
viscoelastic models, while additional dispersion and attenuation effects are produced by
the shear viscosity. Comparisons between the computation times for serial and parallel
implementations show that the nVidia Tesla K20 GPU accelerates these Green’s function
calculations by factors of 45 and 700 for shear wave simulations in elastic and viscoelastic
models, respectively, relative to serial calculations on a desktop computer with an Intel(R)
Core(TM) i7-2600 3.40GHz CPU. When these Green’s functions are evaluated for push
beams with different spatial samplings, the results show that the percent difference is larger
for shear wave particle displacements calculated by push beams with lower f/#. The results
indicate that the percent difference usually decreases monotonically as the spatial sampling
decreases, with the exception of the simulated shear waves in a viscoelastic medium using
f/1.5 and f/2 push beams, where the percent difference for both of these reach a local
minimum when the spatial sampling is λ/2. The results also indicate that effective results
are obtained using push beams sampled at λ/2, which yield percent differences less than 5%.
Measured and simulated shear waves in elastic and viscoelastic shear wave phantoms are
compared, and the results suggest that the Green’s function simulation approach provides
an effective simulation model for shear waves in elastic and viscoelastic shear wave phantoms.

54

CHAPTER 4
TIME-DOMAIN ESTIMATION OF SHEAR ELASTICITY AND SHEAR
VISCOSITY FOR SHEAR WAVE ELASTICITY IMAGING (SWEI) IN
VISCOELASTIC SHEAR WAVE PHANTOMS

4.1

Introduction

Shear wave elasticity imaging (SWEI) detects the mechanical properties of soft tissue
by analyzing the shear wave induced by an acoustic radiation force [42]. In SWEI, one or
more focused or unfocused ultrasound beams are applied to the soft tissue for a short period
of time to generate shear waves. The shear waves are measured by comparing consecutive
in-phase/quadrature (IQ) data obtained with a speckle tracking algorithm [11]. Mechanical
properties are then obtained by analyzing the measured shear waves. SWEI has significant
diagnostic potential to noninvasively characterize tissue mechanical properties, which often
correlate with tissue physiology [42, 61]. For example, studies have shown that there is
a positive correlation between the stiffness of the liver tissue and the stage of the liver
fibrosis [44, 61–65]. Other clinical applications for SWEI include quantitative assessment of
mechanical properties in thyroid [68], kidney [71], and prostate [72].
In SWEI, the shear elasticity is often obtained using the time-of-flight (TOF) of propagating shear waves. In TOF measurements for SWEI, the arrival times of the shear waves
are measured at several successive spatial points located at the same depth but in different
lateral positions. One approach obtains the arrival times by measuring the time to peak
(TTP) [61, 157], which tracks the arrival times of the peaks for each measured shear wave.
Another approach obtains the arrival time delays with cross-correlations [81, 158, 159]. Of
these, the cross-correlation approach is preferred because this approach is more robust and
much less sensitive to noise. The shear wave propagation speed cs is then calculated from the
slope of the linear regression between the lateral locations of the spatial points and the shear
wave arrival times. Another approach that improves the robustness of the linear regression
55

uses random sample consensus (RANSAC), which updates the linear fit by randomly generating trial solutions and iteratively updating the model by eliminating the outliers of the
arrival time delays [160, 161]. The shear elasticity µ is then obtained from µ = ρc2s , where ρ
is the density.
Although the TOF method often yields effective values in elastic (lossless) shear wave
phantoms [61,162–164], this method is based upon some simplifications and approximations
that cause significant bias when applied to mechanical parameter estimates in soft tissue,
which is always lossy (viscoelastic) [79, 116]. In viscoelastic media, the shear wave speed
changes during propagation due to shear wave attenuation and dispersion, so the relationship
µ = ρc2s no longer holds [115, 148]. Thus, the shear elasticity estimates provided by the
TOF method are biased by the shear viscosity. Another problem is that the TOF method
implicitly assumes that the shear wave is a plane wave. However, complicated diffraction
effects are observed in all shear waves induced by a three dimensional (3D) acoustic radiation
force [165]. Also, the TOF method yields different values of the estimated shear elasticity
for different transducers, focal depths, and lateral ranges applied to the same viscoelastic
shear wave phantom [166], which suggests that the diffraction of the push beam and the
diffraction of the shear wave also contribute to the bias observed in shear wave parameter
estimates obtained from the TOF method.
To address the above problems with the TOF method, we have developed a new approach
to estimate the shear elasticity and shear viscosity parameters that considers the effects of
the viscoelasticity of the medium, the diffraction of the acoustic radiation force push, and the
diffraction of the shear wave. After the mathematical theory for shear wave propagation in
viscoelastic media is described, a forward model is provided for the simulation of shear waves
induced by an acoustic radiation force push. This model calculates the time-domain Green’s
function for shear waves propagating in a viscoelastic medium [148] with a GPU-based
parallel algorithm [165]. Then, a new approach for estimating the shear elasticity and shear
viscosity values of the medium is described, where the new approach minimizes the difference

56

between simulated and measured shear wave particle velocities. From measured shear waves
in viscoelastic shear wave phantoms, parameters of the shear elasticity and viscosity are
estimated and quantitative images are generated. Shear elasticity values estimated with the
conventional time-of-flight method are then compared with the values estimated with the
proposed model-based estimation approach. In addition, the results provided by the modelbased estimation approach with different push beams in viscoelastic shear wave phantoms are
shown and compared. The results indicate that this new model-based approach is an effective
method for extracting the values for the shear elasticity and viscosity from viscoelastic shear
wave phantoms.

4.2
4.2.1

Theory
Acoustic radiation force

In shear wave elasticity imaging (SWEI), shear waves in soft tissue are induced by an applied
acoustic radiation force, which is produced by a focused ultrasound beam. The acoustic
radiation force is modeled as [41]
f (r, t) =

2αI(r, t)
,
cp

(4.1)

where f (r, t) is the body force, r denotes a spatial point in Cartesian coordinates, t denotes
time, α is the attenuation coefficient of the compressional wave, cp is the compressional
wave speed, and I(r, t) is the acoustic intensity. The intensity of the ultrasound beam is
proportional to the square of the acoustic pressure, which is defined as I = |p|2 /ρcp . In
SWEI, the acoustic radiation force is applied to the tissue for a short period of time to
produce the shear wave. To account for the temporal component of the acoustic radiation
force, f (r, t) is written as f (r, t) = f (r)w(t), where f (r) is the spatial distribution of the
acoustic radiation force, and w(t) is a window function that describes the time dependence of
the acoustic radiation force. The window function is defined here as w(t) = rect[(t−T /2)/T ],
where rect(·) is the rectangle function, and T is the duration of the push.

57

4.2.2

Shear waves in viscoelastic media

Shear wave propagation in a viscoelastic medium is governed by Navier’s equation combined
with Voigt’s viscosity model [148, 152]




∂
∂
∂
λ + 2µ + (ηp + 2ηs )
∇(∇ · u) + µ + ηs
∇ × (∇ × u) + f = ρ u,
∂t
∂t
∂t

(4.2)

where u is the shear wave particle displacement, f is the external body force, and t is the
time. The bulk modulus and shear modulus are indicated by λ and µ, respectively, the bulk
viscosity and shear viscosity are indicated by ηp and ηs , respectively, and ρ is the density.
This equation models the elastodynamic response to an applied body force in a viscoelastic
medium.
An approximate solution to this equation is given by [148, 165]
gz (r, t) = gzP (r, t) + gzS (r, t),

(4.3)

2 2
 2

2 − (t−r/cp ) cp
1
z
1
3z
1
1
P
2v
t
p
p
+
e
−
1
Ip (r, t),
(4.4)
gz (r, t) =
4πρcp 2πvp t r3
4πρ r2
r3


 2

2 2
s ) cs
1
1
3z
1
1
z 2 1 − (t−r/c
S
2vs t
√
e
+
gz (r, t) =
1− 2
−
1
Is (r, t),
(4.5)
4πρcs 2πvs t
r
4πρ r2
r
r3




 
 

r 2 c2
!
t−
p
r
p
tc2
cp
p
cp t −
 t
νp t  − 2ν
−
2νp t
p −e
 + erf pcp t
e
 p cp  ,
−
erf
Ip (r, t) = √
 2
2νp t
2νp t
2πcp 

(4.6)
√
νs t
Is (r, t) = √
2πcs


2
 − tcs
e 2νs





2
t− crs c2s

−
+
2νs t
−e





t
erf
2



c t
√s
2νs t




 
cs t − crs
 ,
− erf  √
2νs t


(4.7)
where gzP (r, t) and gzS (r, t) represent the compressional and shear terms for the Green’s
function, respectively, and r denotes a point in Cartesian coordinates. The term gzP (r, t)
p
models the propagation of a compressional wave with propagation speed cp = (λ + 2µ)/ρ,
and the shear term gzS (r, t) models the propagation of a shear wave with propagation speed
58

cs =

p
µ/ρ. In Eqs. 4.3-4.7, νp = (ηp + 2ηs )/ρ, νs = ηs /ρ, and erf(·) denotes the error

function.

4.3

Materials and Methods

4.3.1

Shear wave measurements

Shear waves were measured with a Verasonics ultrasound system (Verasonics Inc., Redmond,
WA) and a C5-2v curvilinear transducer array. To induce shear waves at different focal
depths, three different push beams are generated by the transducer array by changing the
number of active elements and the phase of each transducer element. The number of active
elements and the focal depth for each push beam are given in Table 4.1. Each push beam is
applied for 915 µs with a center frequency of 2 MHz.

Push beam 1
Push beam 2
Push beam 3

Active elements
40
60
92

Push beam focal depth
33 mm
45 mm
70 mm

Table 4.1: Parameters for the three push beams.
After the push beam is applied, the Verasonics system is immediately switched to plane
wave imaging mode to measure the shear waves induced by the push beam. In plane wave
imaging mode, multiple frames of in-phase/quadrature (IQ) data within a two-dimensional
(2D) field-of-view (FOV) are obtained for 36.58 ms with a pulse repetition frequency (PRF)
of 1.832 kHz. The 2D FOV ranges from -36.40 mm to 36.40 mm in the lateral direction
and from 0 to 96.59 mm in the axial direction with a spatial sampling of 0.3450 mm in each
direction. The ultrasound center frequency for plane wave imaging is 4.464 MHz. The shear
wave particle velocities are then calculated from the IQ data using one-dimensional (1D)
autocorrelations [11].

59

0

20

20

20

40

40

40

60
80

100
-20
0
20
(a)
x (mm)

z (mm)

0

z (mm)

z (mm)

0

60
80

60
80

100
-20
0
20
(b)
x (mm)

100
-20
0
20
(c)
x (mm)

Figure 4.1: 2D cross-sections of the three push beams simulated in FOCUS. Figures (4.1a),
(4.1b), and (4.1c) show push beams with focal depths of 33 mm, 45 mm, and 70 mm,
respectively.
4.3.2

Shear wave simulations

To simulate the shear waves induced by the push beams, 3D ultrasound fields are calculated, and then the acoustic radiation force is obtained from Eq. 4.1. The 3D pressure fields are calculated in FOCUS, the “Fast Object-Oriented C++ Ultrasound Simulator” (http://www.egr.msu.edu/∼fultras-web/). For these calculations, the fast near field
method [130,131] evaluates the pressure in a plane in front of the transducer array, and then
the angular spectrum approach [132, 133] computes the 3D pressure field. Fig. 4.1 shows
three push beams with push focal depths of 33 mm, 45 mm, and 70 mm (Table 4.1) that are
simulated for a C5-2v curvilinear transducer array. In each simulation, the sound speed for
the compressional wave is 1540 m/s, and the acoustic attenuation coefficient for the compressional wave is 0.5 dB/cm/MHz. The width of each ultrasound transducer element is
0.460 mm, the element height is 13 mm, the kerf between adjacent elements is 0.048 mm,
and the radius of curvature for each array element in the elevation direction is 60 mm. The
radius of curvature of the transducer array is 49.57 mm.
60

The shear wave particle displacements are then calculated by spatially and temporarily
convolving the Green’s function in Eqs. 4.3-4.7 with the acoustic radiation force. This spatiotemporal convolution is given by


Ny Nz
Nt
Nx X
X
X
X



uz (r, t) =
gz r − r0 , t − τ f (r0 )∆x∆y∆z  w(τ )∆τ,
i=1

(4.8)

m=1 n=1 k=1

where uz (r, t) represents the shear wave particle displacement for the spatial coordinate r
at time t. In Eq. 4.8, the triple-summation within the square brackets represents the spatial
convolution, where Nx , Ny , and Nz represents the number of spatial points in each direction, and ∆x, ∆y, and ∆z are the spatial increments in each direction. The summation
outside of the square brackets represents the temporal convolution, where Nt is the number of time samples, and ∆τ is the time sampling interval. After the shear wave particle
displacement is computed, the shear wave particle velocity vz (r, t) is obtained by numerically calculating the temporal derivative of the shear wave particle displacement uz (r, t)
as vz (r, t) = [uz (r, t + ∆t) − uz (r, t)] /∆t, where ∆t is the time sampling interval. The
simulated shear wave particle velocity values are then utilized within the new model-based
approach to estimate the shear elasticity and shear viscosity of the viscoelastic shear wave
phantoms.
To accelerate the simulations that evaluate Eq. 4.8, the calculations are performed on
a Graphics Processing Unit (GPU) [165]. For these GPU calculations, the pre-calculated
acoustic radiation force, which is stored in a data file, is first loaded into the CPU (host)
memory and then copied to the GPU (device) memory. The kernel function is then called
to separate the evaluation of the spatial convolution in Eq. 4.8 (the triple-summation within
the square brackets) into different threads that are distributed to different processors on the
GPU for parallel evaluations. For each thread, the spatial convolution is evaluated at a single
point in space and at a single moment in time. Each thread obtains the push beam from the
global memory, and the calculated spatial convolution result is also stored in global memory.
After all of the threads are finished, the results are copied back into the host memory. The

61

temporal convolution in Eq. 4.8 (the summation outside of the square brackets) is then
evaluated on the CPU to obtain the shear wave particle displacements, and the shear wave
particle velocity is calculated by evaluating the temporal derivative of Eq. 4.8. The simulated
shear wave particle velocities are stored in a data file for subsequent evaluations. These
calculations are implemented through the Compute Unified Device Architecture (CUDA)
programming language on an nVidia (Santa Clara, CA) Kepler K80 GPU.

4.3.3

Phantom experiments

Shear waves are measured in viscoelastic shear wave phantoms provided by Computerized
Imaging Reference Systems (CIRS), Inc (Norfolk, Virginia, US) with a Verasonics ultrasound
system and a C5-2v transducer. Shear waves are measured in three phantoms (E2348-1,
E2348-2, and E2348-3) with different mechanical properties. For each measurement, the
shear wave particle velocities are obtained for all of the observation points in a 2D plane. An
example of the measured shear waves in phantom E2348-1 using a push beam with a focal
depth of 33 mm is shown in Fig. 4.2. Shear wave measurements are performed ten times on
each phantom with all three different push beams.
0

20
40
60
-30
(a)

-15

0
15
x (mm)

30

time = 6ms

0

20

z (mm)

time = 2ms

z (mm)

z (mm)

0

40
60
-30
(b)

-15

0
15
x (mm)

30

time = 10ms

20
40
60
-30
(c)

-15

0
15
x (mm)

30

Figure 4.2: The measured shear wave particle velocity induced by the push beam with a focal
depth of 33 mm in phantom E2348-1. Figures (4.2a), (4.2b), and (4.2c) show the measured
shear wave particle velocity at 2 ms, 6 ms, and 10 ms, respectively.

62

4.3.4

Time-of-flight method for shear elasticity estimation

In the TOF method, shear elasticity is estimated using cross-correlations [178]. To estimate
the shear elasticity value at point (m, n) where m and n are the indices for the x and y
directions, respectively, the cross-correlations of the measured shear wave particle velocities
are calculated for multiple pairs of points which are located at either side of point (m, n) at
the same depth. Those pairs of points are: {(m − 5, n), (m + 1, n)}, {(m − 4, n), (m + 2,
n)}, {(m − 3, n), (m + 3, n)}, {(m − 2, n), (m + 4, n)}, and {(m − 1, n), (m + 5, n)}. The
time delay for the shear wave particle velocities measured at each pair of points are obtained
from the cross-correlations, namely ∆t1 , ∆t2 , ∆t3 , ∆t4 , and ∆t5 . The shear wave speed
is then calculated for each point pair as Vi = 6∆x/∆ti , where i = 1, 2, 3, 4, 5 and ∆x is the
spatial sampling interval in the x direction. The correlation coefficient CCi for each point
pair is used as the weight to combine all of the calculated shear wave speeds to obtain the
P
estimated shear wave speed at point (m, n) as V (m, n) = 5i=1 Vi ·CCi . The shear elasticity
µ at point (m, n) is then obtained from µ = ρV 2 , where ρ is the density.

4.3.5

Model-based approach for shear elasticity and shear viscosity estimation

Shear elasticity and shear viscosity values are estimated from the measured shear wave
particle velocity through the following approach. At each point in space, the reference shear
wave particle velocity is represented by vref (r, t), and the simulated shear wave particle
velocity for a specific combination of µ and ηs values is represented by vµ,ηs (r, t). The
estimated shear elasticity and shear viscosity values are obtained by minimizing the difference
between the reference shear wave particle velocity and the simulated shear wave particle
velocity,
arg min
µ, ηs

kvref (r, t) − vµ,ηs (r, t)k
,
kvref (r, t)k

(4.9)

where k.k represents the L2 norm with respect to the time variable. Eq. 4.9 is then minimized
with a procedure that determines the values of µ and ηs at each point in space. This is

63

achieved by first computing the vµ,ηs (r, t) waveforms with different shear elasticity µ and
shear viscosity ηs values calculated on the GPU as described in section 4.3.2. The value of
the shear elasticity µ ranges from 2 kPa to 12 kPa with an increment of 0.2 kPa, and the
value of the shear viscosity ηs ranges from 0 to 4 Pa·s with an increment of 0.1 Pa·s. Eq. 4.9
is then evaluated through an exhaustive search that minimizes the difference between the
simulated and the measured shear waves. Finally, this procedure is applied to the timedomain waveforms collected at each point in space, and 2D images of µ and ηs are generated
with the results.

4.4
4.4.1

Results
Time of flight parameter estimates evaluated with simulated shear wave
data

With the results obtained from 3D shear wave simulations, shear elasticity values are initially
estimated using the TOF method. Shear wave particle velocities induced by the push beam
shown in Fig. 4.1a), which is focused at a depth of 33 mm, are calculated with shear wave
particle velocities obtained from elastic and viscoelastic simulation models. For the elastic
simulation model, the shear elasticity is 4 kPa and the shear viscosity is 0 Pa·s. For the
viscoelastic simulation model, the shear elasticity is 4 kPa and the shear viscosity is 2 Pa·s.
Both simulation models are homogeneous, linear, and incompressible. The simulated shear
wave particle velocities are evaluated in the xz plane with y = 0, which ranges from -20 mm
to 20 mm in the x direction and from 0 to 80 mm in the z direction with a spatial increment
of 0.2 mm in each direction. After the shear waves are calculated, white noise is added to
the shear wave particle velocity at each point in space, where the average noise intensity is
20% of the peak shear wave particle velocity. After the shear waves are calculated, the shear
wave elasticity is then estimated using the TOF method described in section 2.5.1.
Two shear elasticity images obtained with the TOF method are shown in Fig. 4.3.
Fig. 4.3a) shows the result obtained from the elastic simulation model with µ = 4 kPa

64

Elastic simulation model
Viscoelastic simulation model
10kPa
10kPa
0
0
8kPa
6kPa

40
4kPa

60

80
-20
(a)

6kPa

40
4kPa

60

2kPa

80
-20
(b)

0kPa

0
x (mm)

8kPa

20
z (mm)

z (mm)

20

20

2kPa
0kPa

0
x (mm)

20

Figure 4.3: Shear elasticity images generated by the TOF method using noisy 3D shear wave
particle velocity data obtained from elastic and viscoelastic simulation models. Fig. 4.3a)
shows the estimated shear elasticity for the homogeneous elastic simulation model with µ = 4
kPa and ηs = 0 Pa·s, and Fig. 4.3b) shows the estimated shear elasticity for the homogeneous
viscoelastic simulation model with µ = 4 kPa and ηs = 2 Pa·s.
and ηs = 0. Near the center of the image, the estimated shear elasticity reaches a maximum
value of 141.58 MPa, and in the region between x = −4 mm and x = 4 mm and between
z = 10 mm and z = 65 mm, which is indicated by the red box in Fig. 4.3a), the average elasticity value is 737.09 kPa. Both of these values differ from the actual value by a significant
amount. These results suggest that, in order to obtain realistic values for the shear elasticity
µ, the TOF approach should only be applied some distance from the focal peak. Also, within
the region indicated by the red box in Fig. 4.3a), the TOF method consistently fails to yield
accurate values, so Fig. 4.3a) is thresholded at 10 kPa to emphasize the range of values closer
to 4 kPa. The average value of µ evaluated over the entire image is 89.44 kPa, which is also
quite large relative to the actual value for µ. However, on either side of the push beam a
short distance from the focal peak, the shear elasticity estimates are more reasonable. In
the white boxes in Fig. 4.3 that extend from 6 mm to 18 mm and −6 mm to −18mm in

65

the x direction and from 25 mm to 45 mm in the z direction, the average shear elasticity
value estimated with the TOF approach is 4.07 kPa, which corresponds to an average error
of 1.75% in this region. The standard deviation of the estimated shear elasticity within the
white boxes is 0.042 kPa, so the variation within this region is also relatively small.
Fig. 4.3b) shows the shear elasticity estimates obtained from the viscoelastic simulation
model using the TOF method, where large errors are observed across the entire image. Near
the center of the image, the value of the estimated shear elasticity reaches a maximum value
of 129.34 MPa, and in the region between x = −4 mm and x = 4 mm and between z = 10
mm and z = 60 mm, which is indicated by the red box in Fig. 4.3b), the average elasticity
value obtained with the TOF method is 1.546 MPa. The average value of µ obtained with
the TOF method over the entire image is 485.15 kPa, which differs from the actual value
µ = 4 kPa by a significant amount. The values in Fig. 4.3b) are also thresholded at 10 kPa
to facilitate comparisons with Fig. 4.3a). In the regions indicated by the white boxes shown
in Fig. 4.3b), the estimated shear elasticity is 6.10 kPa, which corresponds to an error of
52.50% in these regions. The standard deviation of the estimated shear elasticity within the
white boxes is 0.70 kPa, which is larger than the standard deviation in the same region when
the TOF method is applied to the simulation results obtained with the elastic model. These
results reveal that, although the TOF method provides accurate shear elasticity estimates
in some regions when applied to the elastic simulation model with µ = 4 kPa and ηs = 0
Pa·s, the TOF method yields significant errors everywhere when applied to the viscoelastic
simulation model with µ = 4 kPa and ηs = 2 Pa·s.

4.4.2

New model-based shear wave parameter estimation approach applied to
simulated shear wave data

Shear elasticity and shear viscosity values are also obtained using the new model-based
estimation approach applied to the same simulated shear wave particle velocities with white
noise added. These are evaluated for the elastic and viscoelastic models described in the

66

previous section. To estimate the shear elasticity and shear viscosity values at each point
in space, the objective function in Eq. 4.9 is minimized as described in section 4.3.5. In
calculations with the objective function in Eq. 4.9, vref (r, t) is the simulated shear wave
particle velocity with noise added, and vµ,ηs (r, t) is the noise-free shear wave particle velocity
calculated on the GPU for a particular combination of µ and ηs values. The shear elasticity
µ ranges from 2 kPa to 12 kPa with an increment of 0.2 kPa, and the shear viscosity ηs
ranges from 0 to 4 Pa·s with an increment of 0.1 Pa·s in the GPU simulations. Fig. 4.4
shows the objective functions for the shear wave particle velocities calculated at (x, z) = (10
mm, 30 mm) with the elastic and viscoelastic simulation models. In Fig. 4.4, the “o” marker
indicates the value obtained when Eq. 4.9 is optimized, and the “x” marker indicates the
estimated value obtained with the TOF method. When applied to the simulated shear wave
data produced by the elastic model, the new model-based approach yields a value of µ =
4.00 kPa for the shear elasticity, which is exactly the same as the elasticity in the simulation
model, while the estimated shear elasticity provided by the TOF method is 4.11 kPa, which
differs slightly from the actual value (by 2.71%). The shear viscosity value estimated with
the new model-based approach is 0 Pa·s, which is also correct since the simulation model
is purely elastic. Fig. 4.4(a) indicates that the objective function for shear waves simulated
with an elastic model has a single global minimum and that the new model-based approach
successfully finds this global minimum, even in the presence of additive noise. The TOF
method yields a value that is close to the global minimum when applied to the simulated
shear wave data produced by the elastic model.
When applied to the noisy simulated shear wave data produced by the viscoelastic model,
the shear elasticity value obtained with the new model-based approach is 4.00 kPa, which is
exactly the same as the elasticity from the simulation model, whereas the estimated shear
elasticity produced by the TOF method is 5.69 kPa, which differ from the actual value by
42.25%. The shear viscosity value obtained from the noisy simulated shear wave data with
the new model-based approach is 1.90 Pa·s, which deviates from the actual value by only

67

Figure 4.4: Objective functions for the shear wave particle velocities calculated at (x, z)
= (10 mm, 30 mm) with the elastic and viscoelastic simulation models. Fig. 4.4a shows
the objective function evaluated for the elastic simulation model, and Fig. 4.4b shows the
objective function evaluated for the viscoelastic simulation model. The “o” markers indicate
the estimates obtained with the new model-based approach, where µ = 4.00 kPa and ηs =
0 in Fig. 4.4a and µ = 4.00 kPa and ηs = 1.90 Pa·s in Fig. 4.4b. The “x” markers indicate
the estimated values obtained with the TOF method, where µ = 4.11 kPa and ηs = 0 in
Fig. 4.4a and µ = 5.69 kPa and ηs = 0 in Fig. 4.4b.
5%. Figs. 4.4(a) and 4.4(b) indicate that the objective function for the viscoelastic shear
wave simulation model also has a single global minimum and that the new model-based
approach successfully identifies the location of the global minimum. However, when applied
to simulated shear wave data obtained from the viscoelastic model, the parameter estimates
produced by the TOF method are located far from the global minimum.
Fig. 4.5 shows the shear elasticity and shear viscosity images produced by the new estimation approach applied to simulated noisy shear wave particle velocities calculated with
elastic and viscoelastic models. Figs. 4.5(a) and 4.5(b) show the shear elasticity image and
shear viscosity image, respectively, for the noisy elastic simulation model with µ = 4 kPa
and ηs = 0 Pa·s. In the center of the image, between x = −4 mm and x = 4 mm and
between z = 15 mm and z = 70 mm, the average elasticity value is 3.58 kPa with a standard
deviation of 1.16 kPa. The average shear viscosity value in this region is 0.69 Pa·s, and the
standard deviation is 1.10 Pa·s. These results indicate that the new model-based approach

68

Shear elasticity
20

4Pa s

20
6kPa

3Pa s

30

40

4kPa

z(mm)

30

50

40

2Pa s

50
2kPa

1Pa s

60

70
-20 -10
(a)
Viscoelatic simulation model
z(mm)

10

8kPa

60
0kPa

0 10 20
x(mm)

10

70
-20 -10
(b)

0Pa s

0 10 20
x(mm)

10

8kPa

20

4Pa s

20
6kPa

3Pa s

30

30

40

4kPa

50

z(mm)

Elastic simulation model
z(mm)

10

Shear viscosity

40

2Pa s

50
2kPa

60

70
-20 -10
(c)

1Pa s

60
0kPa

0 10 20
x(mm)

70
-20 -10
(d)

0Pa s

0 10 20
x(mm)

Figure 4.5: Shear elasticity and shear viscosity images generated with parameters obtained
from the new model-based approach applied to simulated noisy shear wave particle velocities
calculated with the elastic and viscoelastic simulation models. Figs. 4.5(a) and 4.5(b) show
the shear elasticity and shear viscosity images, respectively, for the elastic simulation model
with µ = 4 kPa and ηs = 0 Pa·s, respectively. Figs. 4.5(c) and 4.5(d) show the shear elasticity
and shear viscosity images, respectively, for a viscoelastic simulation model with µ = 4 kPa
and ηs = 2 Pa·s.
also encounters difficulties close to the focal peak. However, on either side of the push beam
a short distance from the focal peak, the shear elasticity estimates are much more accurate.
In the regions that extend from 6 mm to 18 mm and −6 mm to −18mm in the x direction

69

and from 25 mm to 45 mm in the z direction, which are indicated by the white boxes in
Figs. 4.5(a) and 4.5(b), the average shear elasticity value is 3.99 kPa, which corresponds to
an error of 0.25% in this region. The standard deviation of the estimated shear elasticity
within the white boxes is 0.04 kPa, which is significantly smaller than the standard deviation
evaluated in the middle of the image. The average shear viscosity value is 0.0002 Pa·s with
a standard deviation of 0.005 Pa·s, which is very close to the actual shear viscosity value in
the simulation model.
Figs. 4.5(c) and 4.5(d) show the shear elasticity image and shear viscosity image, respectively, for the noisy viscoelastic simulation model with µ = 4 kPa and ηs = 2 Pa·s. In the
center of the image, between x = −4 mm and x = 4 mm and between z = 15 mm and
z = 65 mm, the average shear elasticity value is 3.42 kPa with a standard deviation of 1.39
kPa and the average shear viscosity value is 1.78 kPa with a standard deviation of 1.61 kPa.
However, on either side of the push beam a short distance from the focal peak, the shear
elasticity estimates with the model-based approach are much more accurate. In the regions
that extend from 6 mm to 18 mm and −6 mm to −18mm in the x direction and from 25
mm to 45 mm in the z direction, which are indicated by the white boxes in Figs. 4.5(c) and
4.5(d), the average shear elasticity value is 3.99 kPa, which corresponds to an error of 0.25%
in this region. The standard deviation of the estimated shear elasticity within the white
boxes is 0.05 kPa, which is significantly smaller than the standard deviation evaluated in the
middle of the image. The average shear viscosity value is 1.76 Pa·s, which corresponds to an
error of 12.2% in this region. The standard deviation of the estimated shear viscosity within
the white boxes is 0.14 Pa·s, which is also significantly smaller than the standard deviation
evaluated in the middle of the image.

70

Figure 4.6: The objective function in Eq. 4.9 evaluated for different values of the shear
elasticity and shear viscosity for shear wave particle velocities measured at (x, z) = (10 mm,
30 mm) in phantom E2348-1. The “o” marker indicates the value obtained with the new
model-based estimation approach (µ = 3.60 kPa and ηs = 0.80 Pa·s), and the “x” marker
indicates the value produced by the TOF method (µ = 6.70 kPa and ηs = 0).
4.4.3

TOF and model-based estimates of shear wave parameters for viscoelastic
shear wave phantoms

To compare the new model-based estimation approach with the TOF method using shear
wave measured data from viscoelastic shear wave phantoms, the two methods are applied
to the shear wave particle velocities measured at one point in space generated in phantom
E2348-1 by a push beam with a focal depth of 33 mm. The TOF method is applied to the
measured shear wave particle velocities at (x, z) = (10 mm, 30 mm), and the TOF estimate
of the shear elasticity at this location is 6.70 kPa. The model-based estimation approach is
applied to the measured shear wave particle velocities in the same location, and the shear
elasticity and shear viscosity parameters are estimated by minimizing the difference between
the measured shear wave particle velocity and the simulated shear wave particle velocity
in Eq. 4.9. The estimated shear elasticity value with the model-based approach is 3.60

71

1.5
Measured
Simulated with TOF parameters

Normalized shaer velocity

Normalized shaer velocity

1.5

1

0.5

0

-0.5
(a) 0

5

10

15

time (ms)

Measured
Simulated with estimated parameters

1

0.5

0

-0.5
(b) 0

5

10

15

time (ms)

Figure 4.7: Measured and simulated shear wave particle velocities using shear wave parameter values estimated with the time-of-flight (TOF) method and the new model-based
estimation approach. In Figs. 4.7a and 4.7b, the blue dashed line indicates the measured
shear wave particle velocity induced by the push beam with a focal peak at 33 mm in phantom E2348-1, the black dotted line in Fig. 4.7(a) is the the simulated shear wave particle
velocity with the shear wave parameters estimated with the TOF method (µ = 6.70 kPa and
ηs = 0), and the red solid line in Fig. 4.7(b) is the the simulated shear wave particle velocity
with the shear wave parameters (µ = 3.60 kPa and ηs = 0.80 Pa·s) obtained from the new
model-based estimation approach.
kPa, and the estimated shear viscosity value with the model-based approach is 0.80 Pa·s.
This result indicates that the shear elasticity value obtained with the TOF method differs
by 86.11% from the value estimated with the model-based estimation approach. Fig. 4.6
illustrates the objective function in Eq. 4.9 for different combinations of shear elasticity and
shear viscosity values. In Fig. 4.6, the “o” marker indicates the value obtained with the
model-based approach (µ = 3.60 kPa and ηs = 0.80 Pa·s), and the “x” marker indicates the
value computed with the TOF method (µ = 6.70 kPa and ηs = 0). The objective function
for the shear wave data measured in phantom E2348-1 has a single global minimum, and
the new model-based approach successfully locates this global minimum. As for the noisy
simulated viscoelastic shear wave data, the shear wave parameters estimated with the TOF
method are far from the global minimum.
Fig. 4.7 demonstrates measured and simulated shear wave particle velocities using the

72

shear wave parameter values estimated with the TOF method and the new model-based
estimation approach. The shear wave particle velocity is induced by a push beam with a
focal peak at 33 mm in phantom E2348-1 and is then measured at (x, z) = (10 mm, 30
mm), which is the same (x, z) coordinate that is used in Fig. 4.6. Fig. 4.7(a) shows the
measured shear wave particle velocity and the simulated shear wave particle velocity using
the parameters estimated with the TOF method (µ = 6.70 kPa and ηs = 0). This result shows
that, although the peaks arrive at approximately the same time, the two waveforms differ by
a significant amount. Overall, the percent difference between these two waveforms is 111.5%.
Fig. 4.7(b) shows the measured and simulated shear wave particle velocity and the simulated
shear wave particle velocity using parameters obtained with the new model-based approach
(µ = 3.60 kPa and ηs = 0.80Pa·s). The two waveforms in Fig. 4.7(b) demonstrate much
better alignment, where the percent difference between the two waveforms in Fig. 4.7(b) is
42.1%. This suggests that the new model-based approach provides more effective estimates
of the shear elasticity and shear viscosity than the TOF method when evaluated at this
location in viscoelastic shear wave phantom E2348-1.

4.4.4

Shear elasticity and shear viscosity imaging in viscoelastic shear wave
phantoms

Fig. 4.8 shows the shear elasticity images provided by the TOF method and the shear
elasticity and shear viscosity images provided by the new model-based approach for shear
waves induced by a push beam focused at 45 mm in three different phantoms (E2348-1,
E2348-2, and E-2348-3). In the regions that extend from 6 mm to 18 mm and −6 mm
to −18mm in the x direction and from 40 mm to 55 mm in the z direction, which are
indicated by the white boxes in each of the figures, the new model-based estimation approach
provides relatively homogeneous values for µ and ηs . Figs. 4.8(a)-(c) show the shear wave
parameter images obtained with the TOF method and the new model-based approach for
phantom E2348-1. In Fig. 4.8(a), which shows the shear elasticity image generated by the

73

8

62

16
12

45

8

4

62

0

80
(b) -20

0
20
x (mm)

Model-based
10

(Pa.s)

4

28

3

45

2

4

62

1

0

80
(c) -20

0
20
x(mm)

z(mm)

45

(kPa)

28
z(mm)

12

Model-based
10

0

0
20
x(mm)

16

10

16

10

4

28

12

28

12

28

3

45

8

45

8

45

2

62

4

62

4

62

1

0

80
(e) -20

0

80
(f) -20

80
(d) -20

0
20
x (mm)

0
20
x(mm)

z(mm)

10

z(mm)

0

0
20
x(mm)

10

16

10

16

10

4

28

12

28

12

28

3

45

8

45

8

45

2

62

4

62

4

62

1

0

80
(h) -20

0

80
(i) -20

80
(g) -20

0
20
x (mm)

0
20
x(mm)

z(mm)

Phantom E2348-2
z(mm)

16

28

80
(a) -20

Phantom E2348-3
z(mm)

(kPa)

z(mm)

Phantom E2348-1
z(mm)

TOF method
10

0

0
20
x(mm)

Figure 4.8: Shear wave parameter images generated by the TOF method and the new modelbased approach for phantoms E2348-1, E2348-2, and E2348-3 where the shear wave is induced
by a push beam focused at 45 mm. Figs. 4.8(a)-(c) show the shear elasticity image provided
by the TOF method, the shear elasticity image provided by the new model-based approach,
and the shear viscosity image provided by the new model-based approach evaluated for
phantom E2348-1, respectively. Figs. 4.8(d)-(f) show the shear elasticity image provided by
the TOF method, the shear elasticity image provided by the new model-based approach, and
the shear viscosity image provided by the new model-based approach evaluated for phantom
E2348-2, respectively. Figs. 4.8(g)-(i) show the shear elasticity image provided by the TOF
method, the shear elasticity image provided by the new model-based approach, and the shear
viscosity image provided by the new model-based approach evaluated for phantom E2348-3,
respectively.

74

TOF method, the average shear elasticity value within the white boxes is 5.24 kPa with
a standard deviation of 1.27 kPa. In Fig. 4.8(b), which shows the shear elasticity image
generated by the new model-based approach, the average shear elasticity value within the
white boxes is 3.72 kPa with a standard deviation of 0.23 kPa. These results show that the
shear elasticity value estimated with the TOF method is significantly larger than the value
estimated with the new model-based approach and that the standard deviation is also larger
with the TOF method. In Fig. 4.8(c), which shows the shear viscosity image generated by
the new model-based approach, the average shear viscosity value within the white boxes is
1.04 Pa·s with a standard deviation of 0.13 Pa·s.
Figs. 4.8(d)-(f) show the shear wave parameter images for phantom E2348-2.

In

Fig. 4.8(d), which shows the shear elasticity image generated by the TOF method, the
average shear elasticity value within the white boxes is 9.15 kPa with a standard deviation
of 3.06 kPa. In Fig. 4.8(e), which shows the shear elasticity image generated by the new
model-based approach, the average shear elasticity value within the white boxes is 6.31 kPa
with a standard deviation of 0.38 kPa. These results show that the shear elasticity value
estimated with the TOF method is significantly larger than the value estimated with the
new model-based approach, and the standard deviation is also larger with the TOF method.
In Fig. 4.8(f), which shows the shear viscosity image generated by the new model-based approach, the average shear viscosity value within the white boxes is 1.51 Pa·s with a standard
deviation of 0.27 Pa·s.
Figs. 4.8(g)-(i) show the shear wave parameter images for phantom E2348-3.

In

Fig. 4.8(g), which shows the shear elasticity image generated by the TOF method, the
average shear elasticity value within the white boxes is 12.02 kPa with a standard deviation
of 3.30 kPa. In Fig. 4.8(h), which shows the shear elasticity image generated by the new
model-based approach, the average shear elasticity value within the white boxes is 8.14 kPa
with a standard deviation of 0.47 kPa. These results show that the shear elasticity value
estimated with the TOF method is significantly larger than the value estimated with the

75

new model-based approach, and the standard deviation is also larger with the TOF method.
In Fig. 4.8(i), which shows the shear viscosity image generated by the new model-based approach, the average shear viscosity value within the white boxes is 1.87 Pa·s with a standard
deviation of 0.35 Pa·s.

4.4.5

Shear elasticity and shear viscosity imaging with different push beams

Fig. 4.9 shows the estimated shear wave parameter images generated by the TOF method
and the new model-based approach evaluated in phantom E2348-2 with three different push
beams. Figs. 4.9(a)-(d) show the estimated shear wave parameter images for a push beam
focused at 33 mm. In the regions that extend from 6 mm to 18 mm and −6 mm to −18mm
in the x direction and from 25 mm to 40 mm in the z direction, which are indicated by the
white boxes, the model-based estimation approach provides relatively homogeneous values
for µ and ηs . For the region outside of the white boxes, greater variation is observed in µ
and ηs . This suggests that the model-based estimation approach should only be applied on
either side of the focal peak some distance away from the push beam. In Fig. 4.9(a), which
shows the shear elasticity image generated by the TOF method for a push beam focused
at z = 33 mm, the average shear elasticity value within the white boxes is 9.41 kPa with
a standard deviation of 3.04 kPa. In Fig. 4.9(b), which shows the shear elasticity image
generated by the new model-based approach for a push beam focused at z = 33 mm, the
average shear elasticity value within the white boxes is 6.14 kPa with a standard deviation
of 0.43 kPa. These results show that shear elasticity value estimated with the TOF method
is significantly larger than the value estimated with the new model-based approach, and the
standard deviation is also larger with the TOF method. In Fig. 4.9(c), which shows the
shear viscosity image generated by the new model-based approach for a push beam focused
at z = 33 mm, the average shear viscosity value within the white boxes is 1.45 Pa·s with a
standard deviation of 1.92 Pa·s. In Fig. 4.9(d) shows an image that describes the percent
difference between the measured shear wave particle velocity and the simulated shear wave

76

45

6

62

(kPa)

12

4

% difference
10

200%

28

3

28

150%

45

2

45

100%

1

62

50%

0

80
(d) -20

Model-based
10

45

6

3

62

3

62

0

80
(b) -20

0

80
(c) -20

0
20
x (mm)

0
20
x(mm)

z(mm)

9

z(mm)

28

(Pa.s)

0
20
x(mm)

z (mm)

9

Model-based
10

0
20
x (mm)

0

10

12

10

4

10

200%

28

9

28

9

28

3

28

150%

45

6

45

6

45

2

45

100%

62

3

62

3

62

1

62

50%

0

80
(f) -20

0

80
(g) -20

0

80
(h) -20

80
(e) -20

0
20
x (mm)

0
20
x(mm)

0
20
x(mm)

z (mm)

12

z(mm)

10

z(mm)

0

0
20
x (mm)

12

10

12

10

4

10

200%

28

9

28

9

28

3

28

150%

45

6

45

6

45

2

45

100%

62

3

62

3

62

1

62

50%

0

80
(j) -20

0

80
(k) -20

0

80
(l) -20

80
(i) -20

0
20
x (mm)

0
20
x(mm)

0
20
x(mm)

z (mm)

10

z(mm)

45 mm push beam
z(mm)

12

28

80
(a) -20

70 mm push beam
z(mm)

(kPa)

z(mm)

33 mm push beam
z(mm)

TOF method
10

0

0
20
x (mm)

Figure 4.9: Shear wave parameter images generated by the TOF method and the new modelbased approach for phantom E2348-2. Figs. 4.9(a), 4.9(e), and 4.9(i) show the shear elasticity
images generated by the TOF method with push beams focused at 33 mm, 45 mm, and 70
mm, respectively. Figs. 4.9(b), 4.9(f), and 4.9(j) show the shear elasticity images generated
by the new model-based approach with push beams focused at 33 mm, 45 mm, and 70 mm,
respectively. Figs. 4.9(c), 4.9(g), and 4.9(k) show the shear viscosity images generated by
the new model-based approach with push beams focused at 33 mm, 45 mm, and 70 mm,
respectively. Figs. 4.9(d), 4.9(h), and 4.9(l) show the images of the difference between the
measured shear wave particle velocity and the simulated shear wave particle velocity for push
beams focused at 33 mm, 45 mm, and 70 mm, respectively.
particle velocity with the optimized shear wave parameters obtained from the model-based
approach. The average percent difference within the white boxes is 38.96% with a standard
deviation of 6.44%.
Figs. 4.9(e)-(h) show the shear wave parameter images for a push beam focused at 45 mm.
77

In the regions that extend from 6 mm to 18 mm and −6 mm to −18mm in the x direction
and from 40 mm to 55 mm in the z direction, which are indicated by the white boxes in
all the figures in Fig. 4.9(e)-(h), the model-based approach provides relatively homogeneous
values for µ and ηs . For the region outside of the white boxes, the variation in µ and ηs is
much greater. Fig. 4.9(e) shows the shear elasticity images generated by the TOF method for
a push beam focused at 45 mm. The average shear elasticity value within the white boxes
is 9.15 kPa with a standard deviation of 3.06 kPa. Fig. 4.9(f) shows the shear elasticity
image generated by the new model-based approach for a push beam focused at 45 mm. The
average shear elasticity value within the white boxes is 6.31 kPa with a standard deviation
of 0.38 kPa. These results show that shear elasticity value estimated with the TOF method
is significantly larger than the value estimated with the new model-based approach, and the
standard deviation is also larger with the TOF method. Fig. 4.9(g) shows the shear viscosity
image generated by the new model-based approach for a push beam focused at 45 mm. The
average shear viscosity value within the white boxes is 1.51 Pa·s with a standard deviation
of 0.27 Pa·s. Fig. 4.9(h) contains and image that describes the percent difference between
the measured shear wave particle velocity and the simulated shear wave particle velocity
with shear wave parameters obtained from the model-based approach. The average percent
difference within the white boxes is 36.90% with a standard deviation of 6.06%.
Figs. 4.9(i)-(l) show the shear wave parameter images for a push beam focused at 70 mm.
In the regions that extend from 6 mm to 18 mm and −6 mm to −18mm in the x direction
and from 60 mm to 78 mm in the z direction, which are indicated by the white boxes in all
the figures in Fig. 4.9(i)-(l), the new model-based estimation approach provides relatively
homogeneous values for µ and ηs . For the region outside of the white boxes, the variation in
µ and ηs is much greater. Fig. 4.9(i) shows the shear elasticity images generated by the TOF
method for a push beam focused at 70 mm. The average shear elasticity value within the
white boxes is 9.50 kPa with a standard deviation of 2.81 kPa. Fig. 4.9(j) shows the shear
elasticity image generated by the new model-based approach for a push beam focused at 70

78

mm. The average shear elasticity value within the white boxes is 6.30 kPa with a standard
deviation of 0.38 kPa. These results show that shear elasticity value estimated with the TOF
method is significantly larger than the value estimated with the new model-based approach,
and the standard deviation is also larger with the TOF method. Fig. 4.9(k) shows the shear
viscosity image generated by the new model-based approach for a push beam focused at 70
mm. The average shear viscosity value within the white boxes is 1.49 Pa·s with a standard
deviation of 0.25 Pa·s. Fig. 4.9(l) shows the image for the percent difference between the
measured shear wave particle velocity and the simulated shear wave particle velocity with
the shear wave parameters obtained from the model-based approach. The average percent
difference within the white boxes is 37.95% with a standard deviation of 6.13%.

4.5
4.5.1

Discussion
Comparisons between the TOF method and the new model-based approach: simulated 3D shear wave data

Figs. 4.3, 4.4, and 4.5 show the shear wave parameters estimated with the TOF method and
the new model-based approach. For all of the estimated shear elasticity and shear viscosity
images in Figs. 4.3 and 4.5, the values at the center of the image show significant bias relative
to the actual values in the simulation. This is expected in the center of the image, where the
shear wave is still very close to the peak of the push beam. However, at some distance from
the peak of the push beam (|x| from 6mm to 18mm and z from 25 mm to 45 mm), more
reasonable estimates are obtained. For shear waves generated by the elastic simulation model,
the average shear elasticity value inside of the white boxes provided by the TOF method
shown in Fig. 4.3(a) is 4.07 kPa with a standard deviation of 0.042 kPa. This corresponds
to an average error of only 1.75%. The standard deviation is also small. The average shear
elasticity value estimated with the new model-based approach in Fig. 4.5(a) is 3.99 kPa with
a standard deviation of 0.04 kPa. This corresponds to an average error of only 0.25% with
a small standard deviation. Thus, for the shear waves generated by the elastic simulation

79

model, both methods provide accurate estimates for the shear elasticity. However, for the
shear waves produced by the viscoelastic simulation model, the average shear elasticity value
estimated with the TOF method in Fig. 4.3(b) is 6.10 kPa with a standard deviation of 0.70
kPa. This corresponds to an average error of 52.50% with a higher standard deviation. The
average shear elasticity value estimated with the new model-based approach in Fig. 4.5(c) is
3.99 kPa with a standard deviation of 0.05 kPa. This corresponds to an average error of only
0.25% with a small standard deviation. Thus, when applied to shear waves generated by
the viscoelastic simulation model, the new model-based approach provides an accurate shear
elasticity estimate, but the estimated provided by the TOF method contains significant bias.
In addition, the average shear viscosity value estimated with the new model-based approach
in Fig. 4.5(c) is 1.76 Pa·s, which corresponds to an error of 12.22%. Although the shear
viscosity estimate provided by the new model-based approach is not as accurate as the shear
elasticity estimate, the error and the standard deviation (0.14 Pa·s) are still relatively small.
Fig. 4.4 contains another example that compares the performance of the TOF method and
the model-based approach applied to elastic and viscoelastic simulation results. Fig. 4.4(a)
shows that, for shear waves simulated with the elastic model, both the TOF method and the
new model-based approach provide accurate shear elasticity estimates that are close to the
global minimum of the objective function. However, Fig. 4.4(b) shows that, when applied to
shear waves generated by the viscoelastic simulation model, the new model-based estimation
approach provides accurate shear elasticity and shear viscosity values, but the shear elasticity
estimate provided by the TOF method contains significant bias.

4.5.2

Comparisons between the TOF method and the new model-based approach: viscoelastic shear wave phantoms

Figs. 4.6 and 4.7 contain shear wave parameter estimation results for the TOF method and
the new model-based approach applied to a viscoelastic shear wave phantom at one point
in space. Fig. 4.6 shows that the new model-based approach yields effective results, but the

80

shear elasticity estimate provided by the TOF method contains significant bias. Fig. 4.7
compares the measured shear wave particle velocities to the simulated shear wave particle
velocities obtained from the shear elasticity and shear viscosity values estimated with the
TOF method and the new model-based approach. Fig. 4.7(a) reveals that the simulated
shear wave particle velocity using the µ and ηs values estimated with the TOF method
demonstrates poor agreement with the measured shear wave particle velocity, as evidenced
by the large percent difference between the two waveforms (111.5%). Fig. 4.7(b) shows that
the simulated shear wave particle velocity using the µ and ηs values estimated with the
new model-based approach demonstrates much better alignment with the measured shear
wave particle velocity, which is supported by the significantly smaller value of the percent
difference between the two waveforms (42.1%). Also, the above results indicate that the
new model-based approach consistently outperforms the TOF method, since the µ and ηs
values estimated with the new model-based approach achieve much better alignment with
the measured shear wave particle velocity.

4.5.3

Shear wave parameter estimation in different phantoms

Fig. 4.8 contains images of shear wave parameters generated by the TOF method and the
new model-based approach in a viscoelastic shear wave phantom. For phantom E2348-1, the
average µ value estimated with the TOF method is 5.24 kPa with a standard deviation of 1.27
kPa, and the average µ value estimated with the new model-based approach is 3.72 kPa with
a standard deviation of 0.23 kPa. This indicates that the estimated shear elasticity provided
by the TOF method differs by 40.9% relative to the estimated shear elasticity provided by
the new model-based approach. For phantom E2348-2, the average µ value estimated with
the TOF method is 9.15 kPa with a standard deviation of 3.06 kPa, and the average µ value
estimated with the new model-based approach is 6.31 kPa with a standard deviation of 0.38
kPa. This reveals that the estimated shear elasticity provided by the TOF method differs by
45.0% relative to the estimated shear elasticity provided by the new model-based approach.

81

25

Model-based approach - 33 mm push beam
Model-based approach - 45 mm push beam

Shear Elasticity (kPa)

20

Model-based approach - 70 mm push beam
TOF method - 33 mm push beam

15

TOF method - 45 mm push beam
TOF method - 70 mm push beam

10

5

0

E2348-1

E2348-2

E2348-3

Figure 4.10: Comparison of the shear elasticity values estimated in three CIRS viscoelastic
shear wave phantoms with the new model-based approach and the time-of-flight method.
The three groups of bars represent the results for three different phantoms. For each group
of bars, the first three bars indicate the result obtained with the new model-based approach,
and the last three bars represent the result given by the time-of-flight method. The error bars represent the standard deviation of the estimated elasticity calculated from ten
measurements.
For phantom E2348-3, the average µ value estimated with the TOF method is 12.02 kPa
with a standard deviation of 3.30 kPa, and the average µ value estimated with the new
model-based approach is 8.14 kPa with a standard deviation of 0.47 kPa. This indicates
that the estimated shear elasticity provided by the TOF method differs by 47.7% relative to
the estimated shear elasticity provided by the new model-based approach. Fig. 4.10 shows
a bar plot that compares the estimated shear elasticity for these three phantoms. Fig. 4.10
indicates that, in all of the phantoms, the shear elasticity estimate provided by the TOF
method is significantly higher than the estimate of the shear elasticity provided by the new
model-based approach and that the standard deviation is also much higher with the TOF
method.
Fig. 4.8 shows several images of the shear elasticity and the shear viscosity in three
phantoms. For phantom E2348-1, the average ηs value estimated with the new model-based

82

3
2.5

33 mm push beam
45 mm push beam

Shear Viscosity (Pa.s)

70 mm push beam
2
1.5

1
0.5
0
E2348-1

E2348-2

E2348-3

Figure 4.11: Estimated viscosity values obtained with the model-based approach in three
CIRS viscoelastic shear wave phantoms. Three different groups of bars describe the results
for three different phantoms. The three bars in each group indicate the results generated by
pushes at different depths. The error bars indicate the standard deviation of the estimated
elasticity from ten measurements.
approach is 1.04 Pa·s with a standard deviation of 0.13 Pa·s. For phantom E2348-2, the
average ηs value estimated with the new model-based approach is 1.51 Pa·s with a standard
deviation of 0.27 Pa·s. For phantom E2348-3, the average ηs value estimated with the new
model-based approach is 1.87 Pa·s with a standard deviation of 0.35 Pa·s. Fig. 4.11 also
contains a bar plot that compares the estimated shear viscosity of the three phantoms.

4.5.4

Shear wave parameter estimation with different push beams

Fig. 4.9 contains images of shear wave parameter estimates with the TOF method and the
new model-based approach for three different push beams evaluated in phantom E2348-2.
For the white boxes in Figs. 4.9(a), 4.9(e), and 4.9(i), the average shear elasticity values
estimated for the TOF method with push beams focused at 33 mm, 45 mm, and 70 mm
are 9.41 kPa, 9.15 kPa, and 9.50 kPa, respectively. For the white boxes in Figs. 4.9(b),
83

4.9(f), and 4.9(j), the average shear elasticity values estimated with the new model-based
approach for push beams focused at 33 mm, 45 mm, and 70 mm are 6.14 kPa, 6.31 kPa,
and 6.30 kPa, respectively. The estimated shear elasticity values for all three phantoms and
all three push beams are compared in Fig. 4.10, which shows that for all of the phantoms,
the variance of the shear elasticity estimate provided by the new model-based approach is
smaller than that of the TOF method for all three push beams. In Figs. 4.9(c), 4.9(g), and
4.9(k), the shear viscosity value estimated in the white boxes with the new model-based
approach for push beams focused at 33 mm, 45 mm, and 70 mm are 1.45 Pa·s, 1.51 Pa·s,
and 1.49 Pa·s, respectively. The estimated shear viscosity values for all three phantoms and
all three push beams are compared in Fig. 4.11, which shows that for different phantoms,
the variance of the shear viscosity estimate provided by the new model-based approach is
larger than the variance of the shear elasticity estimate provided by the TOF method. In
Figs. 4.9(d), 4.9(h), and 4.9(l), the average percent difference between the measured and
simulated shear wave particle velocities using the µ and ηs values estimated with the modelbased approach for push beams focused at 33 mm, 45 mm, and 70 mm are 38.96%, 36.90%,
and 37.95%, respectively. This result indicates that for all three push beams, the new modelbased approach yields effective estimates for the shear elasticity and shear viscosity on either
side of the focal peak away from the push beam.

4.6

Conclusion

A new model-based approach is introduced that provides estimates for the shear elasticity
and shear viscosity in viscoelastic shear wave phantoms. The shear wave particle velocity is
first calculated with different combinations of the shear elasticity and shear viscosity on a
GPU, then a new model-based approach is employed to determine the parameter combination that yields the smallest percent difference between the measured and simulated shear
waves. The shear wave parameter estimates provided by the new model-based approach
are compared with the TOF method in simulation models, and the results show that the

84

new model-based approach provide more accurate result than the TOF method when applied to simulation data obtained from a 3D viscoelastic shear wave model. Shear wave
parameters are then estimated in measured shear wave data collected from three different
phantoms for three different push beams. The results indicate that the shear elasticity value
estimated with the TOF method is significantly higher than the shear elasticity value estimated with the new model-based approach, and the standard deviation is also higher with
the TOF approach. The results also indicate that the shear elasticity values estimated with
the new model-based approach consistently demonstrate much smaller standard deviations
for different push beams relative to those produced by the TOF method.

85

CHAPTER 5
SHEAR ELASTICITY AND SHEAR VISCOSITY ESTIMATION IN LIVER

5.1

Introduction

Liver fibrosis and cirrhosis are wound healing responses to chronic liver injury from a
variety of causes including viral, autoimmune, drug induced, cholestatic, and metabolic diseases [179]. Liver cirrhosis affects many patients. In the United States, an estimated 900,000
patients are suffering from liver cirrhosis, and approximately 30,000 deaths are caused by
liver cirrhosis per year [179]. Quantitative bio-markers are critically needed for more accurate
liver fibrosis diagnostics and staging. Recently, transient ultrasound elastography techniques
have proven to effectively provide a positive correlation between the stiffness of the liver tissue and the stage of the liver fibrosis [44, 61–65, 77, 79]. Noninvasive ultrasound techniques
that determine liver stiffness include shear wave elasticity imaging (SWEI) [44, 61–65] and
supersonic shear imaging (SSI) [77, 79].
In SWEI, one or more focused or unfocused ultrasound beams are applied for a short
period of time to generate tissue deformation. In turn, this deformation produces shear
waves. To measure the shear wave particle displacements, multiple frames of radio-frequency
(RF) or in-phase/quadrature (IQ) data within a two-dimensional (2D) field-of-view (FOV)
are obtained with pulse-echo ultrasound sequences. Then, shear wave particle displacements
are calculated from the RF or IQ data using one-dimensional (1D) autocorrelations between
consecutive frames [11, 42, 60]. In SSI, several push beams focused at different depths are
sequentially applied to the tissue. The focus of the push beams moves downward faster than
the propagation speed of the shear waves, hence these are “supersonic” relative to the shear
wave speed. After the sequence of push beams is applied, shear waves are generated that
propagate away from the push location with a cone-shaped wavefront. After the shear waves
are generated, multiple frames of RF or IQ data within a 2D FOV are measured using an

86

ultrafast compound imaging technique [75, 76], which measures ultrasound signals at a very
high frame rate (up to 6000 frames per second). Then, shear wave particle displacements
are calculated using 1D autocorrelations between consecutive frames [11].
In SWEI and SSI, the shear elasticity value, which represents the stiffness of the tissue,
is usually obtained using the time-of-flight (TOF) method. The TOF method calculates the
differences in the arrival times of the shear waves at several spatial points located at the
same depth but in different lateral positions to determine the shear wave propagation speed.
The differences in the arrival times are determined from the time to peak (TTP), which
identifies the peak of each measured shear wave [61, 157], or from cross-correlations that
calculate the arrival time differences [81, 158, 159]. The shear wave propagation speed cs is
then estimated from the slope of the linear fit between the lateral positions of the observation
points and the differences in the arrival times of the shear waves measured at these points.
The shear elasticity µ is then obtained from µ = ρc2s , where ρ is the density. The TOF
method provides accurate shear elasticity estimates in purely elastic media [61, 162–164].
However, because the TOF method neglects the viscosity, this method yields significant bias
when applied to estimates of µ in soft tissue, which is viscoelastic [79,116]. Furthermore, the
TOF method assumes that the shear wave is a plane wave. However, complicated diffraction
patterns are observed in shear waves induced by a three-dimensional (3D) acoustic radiation
force [165], and the TOF method provides different values of the estimated shear elasticity
for different transducers, focal depths, and lateral ranges as a consequence of shear wave
propagation [166]. Thus, these effects are also important when estimating the mechanical
properties of soft tissue.
Other methods have also been developed for shear wave parameter estimation. One of
these is attenuation measuring ultrasound shearwave elastography (AMUSE) [167]. With
the AMUSE method, a k-space representation of the shear wave is calculated at several
successive spatial points located at the same depth but in different lateral positions using
the two-dimensional fast Fourier transform (2D-FFT). The k-space representation of the

87

shear wave data is plotted as function of the angular frequency (ω) and the wavenumber in
the x-direction (kx ). Then, the shear wave propagation speed cs is obtained at each frequency
by determining the location of the peak at each frequency in k-space and dividing the value
of ω by the value of kx . The attenuation αs is obtained by measuring the full-width-at-halfmaximum (FWHM) in the kx direction for each angular frequency ω. The AMUSE method
provides a model-free estimate of shear wave speed and attenuation. However, this method
is only effective in the low-frequency range, and diffraction effects are not considered.
The above problems are addressed with a new model-based approach that estimates
the shear elasticity and shear viscosity of soft tissue. In contract to the frequency-domain
analysis with AMUSE, the new approach estimates tissue mechanical parameters in the
time-domain. In this new approach, shear wave particle velocities are first simulated using
an approximate viscoelastic Green’s function model on the GPU [165]. The shear elasticity
and shear viscosity values are then estimated using a model-based approach that minimizes
the difference between measured and simulated shear wave particle velocities. Shear elasticity and shear viscosity images are generated by evaluating this approach for all of the
spatial points in a two-dimensional (2D) field-of-view (FOV). The new approach is applied
to measured shear waves data from in vivo human livers, and the estimation results are
compared with the TOF method. The results indicate that this new model-based approach
successfully generates shear elasticity and shear viscosity images with the measured shear
wave data from in vivo human livers.

5.2
5.2.1

Theory
Acoustic radiation force

In transient ultrasound elastography, shear waves in biological tissue are induced by an
applied acoustic radiation force (push beam). The acoustic radiation force is generated by

88

a focused ultrasound beam and is modeled as
f (r, t) = f (r)w(t) =

2αI(r)
w(t),
cp

(5.1)

where α is the acoustic attenuation coefficient for the compressional wave, cp is the sound
speed for the compressional wave, r is the spatial coordinate, t is the time, f (r, t) is the
spatio-temporal acoustic radiation force at location r and time t, I(r) is the time-average
ultrasound intensity at location r, and w(t) is the temporal window that describes the time
duration of the applied acoustic radiation force. In clinical ultrasound elastography, the
window function is described by w(t) = rect[(t − T /2)/T ], where rect(·) is the rectangle
function, and T is the push duration.

5.2.2

Navier’s equation

Shear wave propagation in elastic, linear, and isotropic biological tissue is modeled by
Navier’s equation [152]
(λ + 2µ)∇(∇ · u) − µ ∇ × (∇ × u) + f = ρ

∂
u,
∂t

(5.2)

where u represents the particle displacement, f represents the applied acoustic radiation
force, t represents the time, ρ represents the density, λ and µ represent the bulk modulus and
shear modulus, respectively, and ρ represents the density. Eq. 5.2 models the elastodynamic
wave induced by the acoustic radiation force in a purely elastic medium. If the viscoelasticity
of the biological tissue is considered, Voigt’s model is incorporated into to Navier’s equation
[148, 154], yielding




∂
∂
∂
λ + 2µ + (ηp + 2ηs )
∇(∇ · u) + µ + ηs
∇ × (∇ × u) + f = ρ u,
∂t
∂t
∂t
where ηp and ηs represent the bulk viscosity and shear viscosity, respectively.

89

(5.3)

5.2.3

Green’s functions for Navier’s equation

The Green’s function for Navier’s equation evaluated in an elastic medium calculates the
displacement induced by a spatially and temporally impulsive body force. When the acoustic
radiation force push is applied in the z-direction, the z-component of the Green’s function
that describes the resulting displacement in an elastic medium is represented by [152]
gz (r, t) = gzP (r, t) + gzS (r, t) + gzC (r, t), where


1 z2
r
P
,
gz (r, t) =
δ t−
cp
4πρc2p r3


 
1
z2 1
r
S
,
gz (r, t) =
δ t−
1− 2
r
cs
4πρc2s
r

 
 2



1
1
r
3z
r
C
gz (r, t) =
− 1 3t H t −
−H t−
.
4πρ r2
cp
cs
r

(5.4)
(5.5)
(5.6)
(5.7)

In Eqs. 5.4 - 5.7, gzP (r, t), gzS (r, t), and gzC (r, t) represent the compressional, shear, and the
coupling components of the Green’s function for Navier’s equation evaluated in an elastic
medium, respectively, where r denotes a spatial point in Cartesian coordinates, the comp
p
pressional wave speed is cp = (λ + 2µ)/ρ, and the shear wave speed is cs = µ/ρ. In the
above expressions, δ(·) denotes the Dirac delta function, H(·) denotes the Heaviside (unit
step) function, r = |r| is the distance between the source point and the observation point,
and z is the distance between the source point and the observation point in the z-direction.
An approximate Green’s function for Navier’s equation evaluated in a viscoelastic medium

90

is [148]
gz (r, t) = gzP (r, t) + gzS (r, t) + gzC (r, t), where

(5.8)

(t−r/cp )2 c2p

gzP (r, t) =
gzS (r, t) =
gzC (r, t) =

Ip (r, t) =

1
1
z2 −
2vp t
p
e
,
4πρcp 2πvp t r3


2 2
s ) cs
1
1
z 2 1 − (t−r/c
2v
t
s
√
1− 2
e
,
4πρcs 2πvs t
r
r

 2

1
1 
3z
−
1
I
(r,
t)
+
I
(r,
t)
,
p
s
4πρ r2
r3





r 2 c2
2
t−
p
p
tc
cp

νp t  − 2νp
−
2νp t
p −e
e
 + t erf
√

 2
2πcp

(5.9)
(5.10)
(5.11)
cp t
p
2νp t

!


− erf 



r
cp

cp t −
p
2νp t

 
 ,

(5.12)




2



 

 
r
√
t− crs c2s


2
tc
c
t
−
s
s

cs
νs t 
e− 2νs − e− 2νs t  + t erf √cs t
 ,
− erf  √
Is (r, t) = √


2
2νs t
2νs t
2πcs
(5.13)
where νp = (ηp + 2ηs )/ρ, νs = ηs /ρ, and erf(·) denotes the error function.

5.3
5.3.1

Methods
Acoustic radiation force simulations

In SWEI and SSI, shear waves are generated by focused or unfocused ultrasound beams.
The pressure field is calculated first, and then the radiation force is obtained from the
computed pressure field, where the radiation force then provides the input to Eq. 5.2 or
Eq. 5.3 for calculations of the shear wave particle displacements. The 3D acoustic pressure fields are calculated in FOCUS, the “Fast Object-oriented C++ Ultrasound Simulator”
(http://www.egr.msu.edu/∼fultras-web). FOCUS calculates the full 3D acoustic pressure
field with the fast near field method (FNM) [130, 131] combined with the angular spectrum
approach (ASA) [132, 133]. After the 3D acoustic pressure field is computed, the acoustic

91

radiation force is calculated from Eq. 5.1.

5.3.2

Calculating shear wave particle velocities induced by an acoustic radiation
force using Green’s functions

With the models described in Eqs. 5.4-5.6 and Eqs. 5.8-5.13, the shear wave particle displacement induced by an acoustic radiation force is calculated by evaluating the spatio-temporal
convolution between the Green’s function and the push beam as
ZZZ Z
gz (r − r0 , t − τ )f (r0 , τ )dτ dV0 ,
uz (r, t) =
V0

(5.14)

τ

where uz (r, t) is the induced shear wave particle displacement at a spatial point r, gz (r, t)
is the Green’s function for the elastic or viscoelastic model, and f (r0 , t) = f (r0 )w(t) is
the acoustic radiation force. This convolution is evaluated numerically by superposing the
contributions from all points in the 3D push beam with


Ny Nz
Nt
Nx X
X
X
X



uz (r, t) =
gz r − r0 , t − τ f (r0 )∆x∆y∆z  w(τ )∆τ,
i=1

(5.15)

m=1 n=1 k=1

where Nx , Ny , and Nz are the number of spatial points in each direction. Although Eq. 5.15
effectively calculates the shear wave particle displacement induced by the push beam in
elastic or viscoelastic media, this calculation is very time-consuming when the number of
observation points is large, especially when shear waves are evaluated in 3D. In addition,
the computation time increases further if multiple combinations of shear elasticity and shear
viscosity parameters are evaluated. For instance, when calculating the shear wave particle
displacements for N shear elasticity values and M shear viscosity values, Eq. 5.15 is evaluated
N M times. In order to reduce the computational effort, an approximate model for shear
wave calculation in viscoelastic media is needed.

5.3.3

Approximate model for shear wave calculation in viscoelastic media

In the far field where |r| ≈ x, the Green’s function gzV (r, t) for the viscoelastic model
(Eqs. 5.8-5.13) is approximated with a convolution between the Green’s function gzE (r, t)
92

for the elastic model (Eqs. 5.4-5.6) and an attenuation term A(r, t) as
gzV (r, t) ≈ gzE (r, t) ∗t A(r, t),

(5.16)

where A(r, t) is the approximate causal solution to the 1D Stokes wave equation [180], given
by
ρcs
A(x, t) = H(t) √
[F (t, x) + F (t, −x)] ,
4 2πγt




x
(cs t − x)2
.
F (x, t) = 1 +
exp −
cs t
2γc2s t

(5.17)
(5.18)

In Eq. 5.17 and 5.18, H(t) is the Heaviside function, and γ = η/ρc2 .
When calculating the shear wave particle displacement induced by an acoustic radiation
force with the approximate viscoelastic model, Eq. 5.14 is represented by
uVz (r, t)

ZZZ Z
≈
V0

gzE (r − r0 , t − τ )A(r − r0 , t − τ )f (r0 , τ )dτ dV0 .

(5.19)

τ

In the far-field, A(r − r0 , t − τ ) is approximated by A(x0 , t − τ ), so Eq. 5.19 is approximated
by
uVz (r, t)

Z
≈

ZZZ
A(x0 , t − τ )

τ

gzE (r − r0 , t − τ )f (r0 , τ )dτ dV0 .

(5.20)

V0

Since the triple integral in Eq. 5.19 is exactly the same as Eq. 5.14, Eq. 5.19 is represented
by a temporal convolution between the shear wave particle displacement induced by a push
beam with an elastic model and the attenuation term, which is given by
uVz (r, t) ≈ uE
z (r, t) ∗t A(x0 , t).

(5.21)

The expression for uE
z (r, t) described in Eq. 3.18 is calculated with
uE
z

Ny Nz
Nx X
X
X

(r, t) =
hz,T (r − r0 , t) f (r0 )∆x∆y∆z.

(5.22)

m=1 n=1 k=1

In this triple sum, Nx , Ny , and Nz describe the number of spatial points in each direction,
and ∆x, ∆y, and ∆z are the spatial sampling intervals in the x, y, and z directions for
93

the simulated 3D acoustic radiation force. The superscript in uE
z (r, t) indicates that this
quantity is the shear wave particle displacement calculated with an elastic model. In Eq. 5.22,
hz,T (r, t) describes the temporal convolution between the Green’s function for the elastic
model gzE (r, t) and the window function w(t), which is given by
Z
hz,T (r, t) =

S
gzE (r, t − τ )w(τ )dτ = hP
z,T (r, t) + hz,T (r, t),

(5.23)

τ


 2


t − r/cp − T /2
1 z2
1
3z
1
=
−
rect
−1 3
2
3
2
T
4πρ r
4πρcp r
r
"
!

#

 
t − r/cp − T /2
1
r
1 2 r2
t − 2 rect
+
2tT − T 2 H t −
−T
, (5.24)
×
2
T
2
cp
cp





 2
t − r/cs − T /2
1
1
z2 1
1
3z
S
hz,T (r, t) =
rect
−1 3
1− 2
−
2
2
r
T
4πρ r
4πρc
r
r
  2 s







1 r
r
t − r/cs − T /2
1
2
2
×
− (t − T ) rect
2tT − T H t −
−
. (5.25)
2 c2s
T
2
cs

hP
z,T (r, t)

The shear wave particle displacement that is predicted by this approximate viscoelastic
model is calculated with the following steps:
1. Calculate the shear wave particle displacement induced by a push beam with the
elastic model using Eq. 5.22.
2. Obtain the attenuation term from Eq. 5.18.
3. Numerically convolve the shear wave particle displacement calculated with the elastic
model and the attenuation term, which yields the approximate shear wave particle displacement in a viscoelastic medium.
Fig. 5.1 shows an example of a calculation with the approximate viscoelastic model.
Fig. 5.1a) shows a 2D cross section of the 3D push beam, which is the same push beam
shown in Fig. 3.1. This result is computed in FOCUS for a Philips L7-4 linear ultrasound
transducer array focused at 25 mm. The red dot in Fig. 5.1a) indicates the observation point,
which is located at x = 10 mm and z = 22.8 mm. In this simulation, the shear elasticity
is 2.25 kPa, and the shear viscosity is 1.0 Pa·s. Fig. 5.1b) shows the shear wave particle
displacement computed with the elastic model in Eq. 5.22, Fig. 5.1c) shows the attenuation

94

term calculated with Eq. 5.18, and Fig. 5.1d) shows the approximate shear wave particle
displacement that is calculated by convolving the attenuation term with the shear wave

z (mm)

10
20
30
40
-10
0
10
(a) x (mm)

1

Elastic model

0.5

0
0
(b)

10
20
t (ms)

1

Attenuation term

Normalized Displacement

Push Beam

Normalized Displacement

0

Normalized Displacement

particle displacement computed with the elastic model.

0.5

30

0
0
(c)

10
20
t (ms)

30

1

Approximate model

0.5

0
0
(d)

10
20
t (ms)

30

Figure 5.1: The results of the individual calculation steps that produce an approximate
shear wave particle displacement. Fig. 5.1a shows the push beam, and the red dot indicates
the observation point, which is located at (x, z) = (10 mm, 22.8 mm). Fig. 5.1b shows the
shear wave particle displacement simulated with the elastic model Eq. 5.22. Fig. 5.1c shows
the attenuation term calculated with Eq. 5.18. Fig. 5.1d shows the approximate shear wave
particle displacement calculated by convolving the results shown in Fig. 5.1b and Fig. 5.1c.
This approximate viscoelastic model significantly reduces the computation time. For
example, to calculate the shear wave particle displacements for N shear elasticity values and
M shear viscosity values, the spatial convolution in the conventional viscoelastic model in
Eq. 5.15 is evaluated N × M times. However, with the approximate viscoelastic model, the
spatial convolution in Eq. 5.22 is only evaluated N times when calculating the shear wave
particle displacements with the elastic model using different µ values. The approximate
shear wave particle displacements in viscoelastic media are then obtained from temporal
convolutions, which take much less time than spatial convolutions.

5.3.4

Validation of the approximate viscoelastic model

To validate the approximate viscoelastic model, shear wave particle displacements calculated
with the viscoelastic model in Eqs. 5.8-5.13 and with the approximate model in Eq. 5.21 are
compared. Shear wave particle displacements are calculated for the push beam shown in
Fig. 5.1(a) at the eight observation points shown in Fig. 3.1. The depth of each observation
95

point is the same as the depth of the focal peak for the push beam, which is 22.9 mm.
The observation points are equally spaced in the lateral direction at this depth, where the
distance from the peak of the push beam to each observation point in the x-direction ranges
from 5 to 40 mm with a increment of 5 mm in the y = 0 plane. The shear wave particle
displacements calculated with the two models are shown in Fig. 5.2. The percent difference
between the shear wave particle displacements calculated with the two models is evaluated
with
%difference =

kuapprox (r, t) − uvisco (r, t)k
,
kuapprox (r, t)k

(5.26)

where uapprox (r, t) is the shear wave particle displacement calculated with the approximate
model in Eq. 5.21, uvisco (r, t) is the shear wave particle displacement calculated using the
viscoelastic model in Eqs. 5.8-5.13, and k.k represents the L2 norm in time. For the eight
observation points, the percent differences between the shear wave particle displacements
calculated with the two models are 3.07%, 2.97%, 3.10%, 3.28%, 3.14%, 3.35%, 3.40%, and
3.56%. These results indicate that the shear wave particle displacements calculated with the
two models are very similar.
The percent difference is also evaluated for observation points in a 2D plane ranging from
-20 mm to 20 mm in the x-direction and from 0 mm to 40 mm in the z-direction, where the
spatial sampling interval in both direction is 1 mm. A 2D image of the percent difference is
shown in Fig. 5.3, which also contains a contour that indicates where the percent difference
is equal to 5%. For the region close to the push beam (|x| < 2 mm), the percent difference
is much larger because of the far field approximation that is employed in the model. The
percent difference is also much larger for observation points located close to the transducer
(z < 10 mm). However, when |x| > 2 mm, and when the depth of the observation point
is similar to the depth of the peak of the push beam, the percent difference is usually less
than 5%. This suggests that the approximate model is an effective approximation for the
viscoelastic model on either side of the push beam at the depth of the focal peak.

96

Normalized Displacement

1
3D Viscoelastic Model
Approximate Model

0.8
0.6
0.4
0.2
0
-0.2
0

10

20

30

40

50

time (ms)
Figure 5.2: Comparison between the shear wave particle displacements simulated with the 3D
viscoelastic model and approximate model. The blue lines indicate the shear wave particle
displacements calculated with the viscoelastic model, and the red dashed lines represent
shear wave particle displacements computed with the approximate model.
5.3.5

Shear wave simulations with the approximate model using GPU

To accelerate shear wave simulations with Eq. 5.21, these calculations are performed on a
Graphics Processing Unit (GPU) [165]. For each calculation on the GPU, the pre-calculated
acoustic radiation force, which is stored in a data file, is first loaded into the CPU (host)
memory and then copied to the GPU (device) memory. The kernel function is then used
to separate the evaluation of the spatial convolution in Eq. 5.22 into different threads that
are distributed to different processors on the GPU for parallel evaluations. For each thread,
the spatial convolution is evaluated at a single point in space and at a single moment in
time. Each thread obtains the push beam from the global memory, and the calculated
spatial convolution result is also stored in global memory. After all of the threads are
finished, the results are copied back into the host memory. The temporal convolution in
97

10%

0

8%

z (mm)

10

6%

20
4%

30

2%

40
-20

0%

-10

0

10

20

x (mm)
Figure 5.3: A 2D image showing the percent difference between shear wave particle displacements calculated with the approximate model in Eq. 5.21 and the viscoelastic model in
Eqs. 5.14. The black contour represents a value of 5% for the percent difference.
Eq. 5.21 is then evaluated on the CPU to obtain the shear wave particle displacements,
and the shear wave particle velocity is calculated by evaluating the temporal derivative of
uVz (r, t) in Eq. 5.21. The simulated shear wave particle velocities are stored in a data file for
subsequent evaluations. These calculations are implemented through the Compute Unified
Device Architecture (CUDA) programming language on an nVidia (Santa Clara, CA) Kepler
K80 GPU.

5.3.6

Tissue mechanical parameter estimation using a genetic algorithm

The above model provides the capability to rapidly simulate shear wave particle velocities
induced by an acoustic radiation force, which then enables the estimation of shear wave
parameters from the measured shear wave particle velocities. The shear elasticity and shear
viscosity values are estimated by minimizing the difference between the simulated and measured shear waves described by
arg min
µ, ηs

kvmeasured (r, t) − vµ,ηs (r, t)k
,
kvmeasured (r, t)k
98

(5.27)

where vmeasured (r, t) is the measured shear wave particle velocity, and vµ,ηs (r, t) is the shear
wave particle velocity simulated with the approximate model in Eq. 5.21.
In chapter 4, this minimization problem is solved by exhaustively searching the 2parameter space to obtain the global minimum. In this chapter, the global minimum of
the objective function in Eq. 5.27 is determined with a genetic algorithm. This is achieved
through the following steps:
1. Calculate the acoustic radiation force as described in section 5.3.1.
2. Calculate the shear wave particle velocity induced by the acoustic radiation force with
the elastic model using Eqs. 5.4-5.6 as described in section 5.3.5. Shear wave particle velocities are calculated with shear elasticities ranging from 0.5 kPa to 4 kPa with an increment
of 0.1 kPa on GPU and are then stored.
3. Use a genetic algorithm to find the global minimum of the objective function in
Eq. 5.27. First, the pre-calculated shear particle velocities simulated with the elastic model
in step 2 are loaded into memory. Then, the minimization problem in Eq. 5.27 is solved
using a genetic algorithm to find the µ and ηs values that minimize the difference between
the simulated and the measured shear wave particle velocities. In Eq. 5.27, vµ,ηs (r, t) is
calculated using the approximate viscoelastic model in Eq. 5.21, which evaluates a temporal
convolution between the pre-calculated shear particle velocities simulated with the elastic
model and the attenuation term. For the optimization, the shear elasticities range from 0.5
kPa to 4 kPa with an increment of 0.1 kPa, the shear viscosities range continuously from 0.1
Pa·s to 4 Pa·s, the population size for the genetic algorithm is 30, and the stopping criteria is
0.1%, which means that if the difference between the best µ and ηs values for the consecutive
generations are less than 0.1%, the optimization will stop. The genetic algorithm is applied
implemented through the Global Optimization Toolbox in MATLAB (Mathworks, MA).

99

0

20

20

40

z (mm)

z (mm)

0

40

60

60

80

80

-10 0 10
(a)
x (mm)

-10 0 10
(b)
x (mm)

Figure 5.4: Two push beams that induce shear waves in human liver. Fig. 5.4a shows the
push beam focused at 35 mm, and Fig. 5.4b shows the push beam focused at 50 mm.
5.3.7

Validation with In Vivo human data

Human liver data was collected at Mayo Clinic under the approval of the Institutional Review
Board at Mayo Clinic in accordance with the Health Insurance Portability and Accountability
Act. Shear wave measurements were performed on healthy subjects with written consent.
The Verasonics Vantage ultrasound system was used in these experiments, and a Philips
P4-2 linear transducer array was used to generate and measure the shear waves, where the
array contains 64 transducer elements with an approximate elevation focus of 70 mm. The
height of each transducer element is 13 mm, the width is 315 µm, and the kerf between
adjacent transducer elements is 20 µm. During each measurement, the ultrasound probe
was inserted into one of the intercostal spaces above the liver, and a homogeneous region of
the liver was located with B-mode ultrasound imaging. Then, shear waves were induced and
measured within that region. Figure 5.4 shows the two push beams that are used to induce
the shear waves. The first push beam in Figure 5.4(a) is electronically focused at 35 mm and

100

is excited at 2.8 MHz. The second push beam in Figure 5.4(b) is electronically focused at 50
mm and is excited at 2.7 MHz. During the experiment, both push beams are excited three
times consecutively for a duration of 267 µs with an interval of 15 µs. Immediately after the
push sequences are completed, the Verasonics ultrasound system is switched to pulse-echo
imaging mode, in which focused ultrasound pulses with a center frequency of 2.78 MHz are
transmitted and the pulse-echo IQ data for the 2D lateral plane are quickly collected. For
the push beam focused at 35 mm, the IQ data is measured for 15.7 ms with a repetition
frequency of 4.78 kHz. For the push beam focused at 50 mm, the IQ data is measured for
17.6 ms with a repetition frequency of 4.27 kHz. Shear wave displacement are then calculated
from the IQ data using the 1D autocorrelation method [11]. For both push beams, shear
wave particle velocities are measured in two different locations: the 7th intercostal space and
the 8th intercostal space. At each location, the measurement is repeated five times. Thus,
there are 2 (push beams) × 2 (positions) × 5 (repetition times) = 20 datasets in total .

5.4
5.4.1

Result
Shear wave particle velocities measured in liver

Fig. 5.5 shows the shear wave particle velocities measured in human liver. Figs. 5.5(a)-5.5(c)
show the shear wave particle velocities induced by the push beam focused at 35 mm at 2
ms, 4 ms, and 6 ms, respectively. All three figures are normalized to the maximum value in
Fig. 5.5(a). In Fig. 5.5(a), the maximum value is located at (x, z) = (3.35 mm, 31.54 mm),
and the amplitude is normalized to 1. At x = 3.35 mm, the full width at half maximum
(FWHM) along the z-direction is 14.09 mm. In Fig. 5.5(b), the maximum value is located at
(x, z) = (5.72 mm, 30.36 mm), and the normalized amplitude is 0.61. At x = 5.72 mm, the
FWHM along the z-direction is 15.45 mm. In Fig. 5.5(c), the maximum value is located at
(x, z) = (8.08 mm, 28.39 mm), and the normalized amplitude is 0.31. At x = 8.08 mm, the
FWHM along the z-direction is 15.21 mm. Also, the amplitude of the shear waves decreases
as the shear wave propagates due to the combined effects of diffraction and attenuation in

101

time = 2ms

0

0
10

20

20

20

30
40

60
(a) -20 -10 0 10 20
x (mm)
0

z (mm)

10

50

30
40
50

30
40

60
(c) -20 -10 0 10 20
x (mm)
0

10

10

20

20

20

40
50

60
(d) -20 -10 0 10 20
x (mm)

z (mm)

10

30

time = 6ms

50

60
(b) -20 -10 0 10 20
x (mm)
0

z (mm)

z (mm)

time = 4ms

10
z (mm)

z (mm)

0

30
40
50

30
40
50

60
(e) -20 -10 0 10 20
x (mm)

60
(f) -20 -10 0 10 20
x (mm)

Figure 5.5: The shear wave particle velocities measured in human liver. Figs. 5.5(a)-5.5(c)
show the shear wave particle velocities induced by the push beam focused at 35 mm at 2 ms,
4 ms, and 6 ms, respectively, and Figs. 5.5(d)-5.5(f) show the shear wave particle velocities
induced by the push beam focused at 50 mm at 2 ms, 4 ms, and 6 ms, respectively.
viscoelastic soft tissue.
Figs. 5.5(d)-5.5(f) show the shear wave particle velocities induced by the push beam
focused at 35 mm at 2 ms, 4 ms, and 6 ms, respectively. All three figures are normalized
to the maximum value in Fig. 5.5(d). In Figs. 5.5(d), the maximum value is located at
(x, z) = (3.35 mm, 38.24 mm), and the amplitude is normalized to 1. At x = 3.35 mm, the
FWHM along the z-direction is 24.24 mm. In Fig. 5.5(e), the maximum value is located at
(x, z) = (5.32 mm, 38.64 mm), and the normalized amplitude is 0.63. At x = 5.32 mm, the
FWHM along the z-direction is 24.12 mm. In Fig. 5.5(f), the maximum value is located at
102

(x, z) = (8.48 mm, 35.09 mm), and the normalized amplitude is 0.35. At x = 8.48 mm, the
FWHM along the z-direction is 22.12 mm. Again, the amplitude of the shear waves decreases
as the shear wave propagates due to the combined effects of diffraction and attenuation in
viscoelastic soft tissue. Comparing Figs. 5.5(d)-5.5(f) with Figs. 5.5(a)-5.5(c), the results
show that the maximum location of the shear wave particle velocity along the z direction
induced by the push beam focused at 35 mm is significantly shallower than the maximum
location of the shear wave particle velocity induced by the push beam focused at 50 mm.
Also, the FWHM of the shear wave particle velocity along the z direction induced by the
push beam focused at 35 mm is also smaller than the FWHM of the shear wave particle
velocity induced by the push beam focused at 50 mm.

5.4.2

Shear elasticity and shear viscosity estimation with the new model-based
approach at a single point in space

After the shear wave particle velocities are measured in the liver, shear wave parameters are
estimated with the new model-based approach. Fig. 5.6 shows a comparison between the
measured shear wave particle velocity and the simulated shear wave particle velocity using
the shear elasticity and shear viscosity values estimated with the new model-based approach.
The measured shear wave particle velocities shown in Fig. 5.6 are induced by the push beam
focused at 35 mm, where Figs. 5.6(a)-5.6(d) show the results obtained at different observation
points. Fig. 5.6(a) shows the shear wave measured at (x, z) = (5 mm, 30 mm). The shear
wave parameters estimated with the new model-based approach in this location are µ = 1.5
kPa and ηs = 1.21 Pa·s. The percent difference between the measured and simulated shear
waves in Fig. 5.6(a) is 29.50%. Fig. 5.6(b) shows the shear wave measured at (x, z) = (6
mm, 35 mm). The shear wave parameters estimated with the new model-based approach
in this location are µ = 1.4 kPa and ηs = 1.31 Pa·s. The percent difference between the
measured and simulated shear waves in Fig. 5.6(b) is 44.11%. Fig. 5.6(c) shows the shear
wave measured at (x, z) = (3 mm, 30 mm) mm. The shear wave parameters estimated with

103

0.5

0

-0.5
(a) 0

5

10
time (ms)

15

1
Measured
Simulated

0.5

0

-0.5
(b) 0

5

10
time (ms)

15

Normalized shear particle velocity

Measured
Simulated

Normalized shear particle velocity

Normalized shear particle velocity
Normalized shear particle velocity

1

1
Measured
Simulated

0.5

0

-0.5
(b) 0

5

10
time (ms)

15

1
0.5
0
-0.5
-1
Measured
Simulated

-1.5
-2
(d) 0

5

10
time (ms)

15

Figure 5.6: Comparisons between measured shear wave particle velocities and shear wave
particle velocities simulated with the µ and ηs values estimated using the new model-based
approach at different locations in space. The blue solid lines indicate the measured shear
wave particle velocities, and the red dashed lines represent the simulated shear wave particle
velocities. Fig. 5.6(a) shows the measured and simulated shear wave particle velocities at
(x, z) = (5 mm, 30 mm), Fig. 5.6(b) shows the measured and simulated shear wave particle
velocities at (x, z) = (6 mm, 35 mm), Fig. 5.6(c) shows the measured and simulated shear
wave particle velocities at (x, z) = (3 mm, 30 mm), and Fig. 5.6(d) shows the measured and
simulated shear wave particle velocities at (x, z) = (15 mm, 30 mm).
the new model-based approach in this location are µ = 4.0 kPa and ηs = 0.02 Pa·s. The
percent difference between the measured and simulated shear waves in Fig. 5.6(c) is 72.30%.
Fig. 5.6(d) shows the shear wave measured at (x, z) = (15 mm, 30 mm). The shear wave

104

3

30

2

40

1

50
-10 0
10
(a) x (mm)

0

Optimized
20

(kPa)

30

4
3
2

40
50
-10 0
10
(b) x(mm)

1
0

Optimized
20

s

(Pa s)
3

30

2

40

1

50
-10 0
10
(c) x(mm)

0

Waveform difference
200%
20
z (mm)

4

z(mm)

(kPa)

z(mm)

z (mm)

TOF method
20

30

150%
100%

40
50
-10 0
10
(d) x (mm)

50%
0

Figure 5.7: Shear wave parameter images for human liver produced by a push beam focused
at 35 mm. Fig. 5.7(a) shows the shear elasticity image generated by the TOF method.
Fig. 5.7(b) shows the shear elasticity image generated by the new model-based approach.
Fig. 5.7(c) shows the shear viscosity image generated by the new model-based approach.
Fig. 5.7(d) shows the percent difference between the measured shear waves and the simulated
shear waves with the µ and ηs values estimated with the new model-based approach.
parameters estimated with the new model-based approach in this location are µ = 1.0 kPa
and ηs = 0.32 Pa·s. The percent difference between the measured and simulated shear waves
in Fig. 5.6(d) is 138.11%.
Figs. 5.6(a)-5.6(d) suggest that the performance of the new model-based approach differs
for different locations in space. For Figs. 5.6(a) and 5.6(b), the shear wave particle velocities
are collected near the peak of the push beam. In these locations, the difference between the
measured and simulated shear wave particle velocities is relatively small. In Fig. 5.6(c), the
shear wave particle velocity is measured too close to the push beam, and in Fig. 5.6(d), the
shear wave particle velocity is measured too far from the push beam. These results suggest
that, in order to obtain effective estimates of the shear elasticity and shear viscosity, the
shear wave particle velocities should be measured in the vicinity of (but not too close to)
the push beam near the depth of the peak of the push beam.

5.4.3

Shear elasticity and shear viscosity imaging in human liver

The new model-based approach is evaluated at all points in space, and images of the shear
wave parameters µ and ηs are generated. Fig. 5.7 shows the shear wave parameter images
for human liver produced by a push beam focused at 35 mm. Fig. 5.7(a) shows the shear
105

elasticity image generated by the TOF method, Fig. 5.7(b) shows the shear elasticity image
generated by the new model-based approach, Fig. 5.7(c) shows the shear viscosity image
generated by the new model-based approach, and Fig. 5.7(d) shows the percent difference
between the measured shear waves and the simulated shear waves with the µ and ηs values
estimated with the new model-based approach. In Figs. 5.7(b)-5.7(d), a contour is superimposed on each figure to indicate the region where:
1. The percent difference shown in Fig. 5.7(d) is less than 60%,
2. The distance from the observation point to the central axis of the push beam is between
3 mm and 10 mm,
3. The estimated shear elasticity values shown in Fig. 5.7(b) are not equal to the lower
bound (1 kPa) or the upper bound (4 kPa), and
4. The estimated shear viscosity value shown in Fig. 5.7(c) is between 0.2 Pa·s and 2.8 Pa·s.
The resulting contour on the left side of the push beam extends from -6.90 mm to -3.35
mm in the x direction and from 24.44 mm to 39.22 mm in the z direction. The contour
on the right side of the push beam spans from 3.75 mm to 6.90 mm in the x direction and
from 21.29 mm to 42.97 mm in the z direction. The total area within the contour is 103.27
mm2 . In Fig. 5.7(b), the average shear elasticity value estimated with the new model-based
approach within the contour is 1.51 kPa with a standard deviation of 0.24 kPa. In comparison, Fig. 5.7(a) shows the shear elasticity value estimated with the TOF method, and the
average estimated shear elasticity value within the same region is 3.92 kPa with a standard
deviation of 6.89 kPa. These results suggest that the new model-based approach provides
a more robust shear elasticity estimate in human liver than the TOF method. Also, the
value for µ obtained with the new model-based approach is significantly smaller than the
value estimated with the TOF method. This is consistent with the behavior observed in
chapter 4. In Fig. 5.7(c), the average shear viscosity value estimated with the new modelbased approach within the contour is 0.94 kPa with a standard deviation of 0.26 kPa. In
Fig. 5.7(d), the average percent difference between the measured and simulated shear wave

106

4

30

3

40

2

50
60
-10 0
10
(a) x (mm)

Optimized
20

(kPa)

4

Optimized
20

s

(Pa s)

30

3

40

2

1

50

1

50

0

60
-10 0
10
(b) x(mm)

0

60
-10 0
10
(c) x(mm)

z(mm)

30

3
2

40
1
0

Waveform difference
200%
20
z (mm)

(kPa)

z(mm)

z (mm)

TOF method
20

30

150%

40

100%

50

50%

60
-10 0
10
(d) x (mm)

0

Figure 5.8: Shear wave parameter images for human liver produced by a push beam focused
at 50 mm. Fig. 5.8(a) shows the shear elasticity image generated by the TOF method.
Fig. 5.8(b) shows the shear elasticity image generated by the new model-based approach.
Fig. 5.8(c) shows the shear viscosity image generated by the new model-based approach.
Fig. 5.8(d) shows the percent difference between the measured shear waves and the simulated
shear waves with the µ and ηs values estimated with the new model-based approach.
particle velocity with the µ and ηs values estimated with the model-based approach is 49.3%
within the contour with a standard deviation of 6.14%.
Fig. 5.8 contains shear wave parameter images generated from human liver data for a push
beam focused at 50 mm. Fig. 5.8(a) shows the shear elasticity image generated by the TOF
method, Fig. 5.8(b) shows the shear elasticity image generated by the new model-based
approach, Fig. 5.8(c) shows the shear viscosity image generated by the new model-based
approach, and Fig. 5.8(d) shows the percent difference between the measured shear waves
and the simulated shear waves using the µ and ηs values estimated with the new model-based
approach. In Figs. 5.8(b)-5.8(d), a contour is again superimposed upon each figure using
the criteria described in the previous paragraph. The contour on the left side of the push
beam extends from -9.26 mm to -3.76 mm in the x direction and from 27.20 mm to 59.14
mm in the z direction. On the right side of the push beam, the larger contour expands from
4.53 mm to 8.87 mm in the x direction and from 27.60 mm to 46.52 mm in the z direction.
Another smaller contour on the right side of the push beam extends from 3.35 mm to 5.32
mm in the x direction and from 51.65 mm to 59.14 mm in the z direction. The total area
within the three contours is 182.84 mm2 . In Fig. 5.8(b), the average shear elasticity value

107

Measurement
Measurement
Measurement
Measurement
Measurement

1
2
3
4
5

Push beam focused at 35 mm
7th rib-cage gap 8th rib-cage gap
1.50 ± 0.25
2.31 ± 0.44
1.54 ± 0.20
2.10 ± 0.32
1.55 ± 0.29
1.99 ± 0.50
1.63 ± 0.22
1.98 ± 0.38
1.68 ± 0.23
1.88 ± 0.29

Push beam focused at 50 mm
7th rib-cage gap 8th rib-cage gap
2.11 ± 0.39
1.82 ± 0.37
2.23 ± 0.33
1.85 ± 0.34
2.28 ± 0.35
1.90 ± 0.29
2.41 ± 0.44
1.91 ± 0.32
2.27 ± 0.38
2.27 ± 0.49

Table 5.1: The average and the standard deviation of the estimated shear elasticity values
provided by the new model-based approach for 20 different measurements. The units are
kPa for each entry in the table.
estimated with the new model-based approach within these contours is 1.85 kPa with a
standard deviation of 0.34 kPa. In contrast, the shear elasticity value estimated with the
TOF method in Fig. 5.8(a) within the same region is 2.46 kPa with a standard deviation
of 0.84 kPa. Again, µ values obtained with the new model-based approach are significantly
smaller than the µ values estimated with the TOF method. This is also consistent with the
behavior observed in chapter 4. In Fig. 5.8(c), the average shear viscosity value estimated
with the new model-based approach within the three contours is 0.94 kPa with a standard
deviation of 0.28 kPa. In Fig. 5.8(d), the average percent difference between the measured
and simulated shear wave particle velocity within the three contours using the µ and ηs
values estimated with the new model-based approach is 52.0% with a standard deviation of
7.97%.

5.4.4

Shear wave parameter estimates in liver

Shear waves are induced and measured in the liver with the transducer placed in 2 different
rib cage gaps with 2 different push beams applied, and for each location and push beam,
5 measurements are repeated. Thus, in total, 20 measurements are performed. In each
measurement, shear wave parameters are collected using the new model-based approach and
the TOF method, and the shear wave parameter values for each measurement are determined
within the contours defined with the criteria described in section 5.4.3. Then, the average

108

and standard deviation are calculated within these contours for each measurement. Table 5.1
shows the shear elasticity values estimated with the new model-based approach for all of these
measurements. For the five shear wave measurements in the 7th rib cage gap with a push
beam focused at 35 mm, the average estimated shear elasticity values are 1.50 kPa, 1.54 kPa,
1.55 kPa, 1.63 kPa, and 1.68 kPa. The standard deviations of the estimated shear elasticity
values within the contours for these five measurements cases are 0.25 kPa, 0.20 kPa, 0.29 kPa,
0.22 kPa, and 0.23 kPa. For the five shear wave measurements in the 8th rib cage gap with
a push beam focused at 35 mm, the average estimated shear elasticity values are 2.31 kPa,
2.10 kPa, 1.99 kPa, 1.98 kPa, and 1.88 kPa. The standard deviation of the estimated shear
elasticity values within the contours for these five measurements are 0.44 kPa, 0.32 kPa, 0.50
kPa, 0.38 kPa, and 0.29 kPa. For the five shear wave measurements in the 7th rib cage gap
with a push beam focused at 50 mm, the average estimated shear elasticity values are 2.11
kPa, 2.23 kPa, 2.28 kPa, 2.41 kPa, and 2.27 kPa. The standard deviation of the estimated
shear elasticity values within the contours for these five measurements are 0.39 kPa, 0.33
kPa, 0.35 kPa, 0.44 kPa, and 0.38 kPa. For the five shear wave measurements in the 8th rib
cage gap with a push beam focused at 50 mm, the average estimated shear elasticity values
are 1.82 kPa, 1.85 kPa, 1.90 kPa, 1.91 kPa, and 2.27 kPa. The standard deviation of the
estimated shear elasticity values within the contours for these five measurements are 0.37
kPa, 0.34 kPa, 0.29 kPa, 0.32 kPa, and 0.49 kPa. All of the above results indicate that the
estimated shear elasticity values estimated with the new model-based approach have a small
variance for each test case within the contours, and the average shear elasticity values also
demonstrate a small amount of variation. However, for different measurement locations and
for different push beams, the differences in the shear elasticity are somewhat larger.
Table 5.2 shows the shear elasticity values for the same measurements estimated with
the TOF method. For the five shear wave measurements in the 7th rib cage gap with a push
beam focused at 35 mm, the average estimated shear elasticity values are 3.92 kPa, 2.33 kPa,
2.33 kPa, 2.33 kPa, and 2.78 kPa. The standard deviations of the estimated shear elasticity

109

Measurement
Measurement
Measurement
Measurement
Measurement

1
2
3
4
5

Push beam focused at 35 mm
7th rib-cage gap 8th rib-cage gap
3.92 ± 6.89
4.04 ± 1.45
2.33 ± 0.84
3.18 ± 0.99
2.33 ± 0.74
3.44 ± 1.26
2.78 ± 1.45
4.51 ± 2.63
3.23 ± 0.95
2.93 ± 1.07

Push beam focused at 50 mm
7th rib-cage gap 8th rib-cage gap
3.57 ± 2.15
2.68 ± 0.81
3.62 ± 1.40
2.46 ± 0.84
3.60 ± 1.41
3.28 ± 1.77
3.80 ± 2.24
2.74 ± 1.33
3.56 ± 1.71
3.56 ± 1.60

Table 5.2: The average and the standard deviations of the shear elasticity values estimated
with the TOF method for 20 different measurements. The units are kPa for each entry in
the table.
values within the contours for these five measurements are 6.89 kPa, 0.84 kPa, 0.74 kPa, 1.45
kPa, and 0.95 kPa. The results show that these standard deviations are much larger than
the standard deviations obtained with the new model-based approach. For the five shear
wave measurements in the 8th rib cage gap with a push beam focused at 35 mm, the average
estimated shear elasticity values are 4.04 kPa, 3.18 kPa, 3.44 kPa, 4.51 kPa, and 2.93 kPa.
The standard deviations of the estimated shear elasticity values within the contours for these
five measurements are 1.45 kPa, 0.99 kPa, 1.26 kPa, 2.63 kPa, and 1.07 kPa. These are also
much larger than the standard deviations obtained with the new model-based approach. For
the five shear wave measurements in the 7th rib cage gap with a push beam focused at 50
mm, the average estimated shear elasticity values are 3.57 kPa, 3.62 kPa, 3.60 kPa, 3.80 kPa,
and 3.56 kPa. The standard deviations of the estimated shear elasticity values within the
contours for these five measurements are 2.15 kPa, 1.40 kPa, 1.41 kPa, 2.24 kPa, and 1.71
kPa. These are once again much larger than the standard deviations obtained with the new
model-based approach. For the five shear wave measurements in the 8th rib cage gap with
a push beam focused at 50 mm, the average estimated shear elasticity values are 2.68 kPa,
2.46 kPa, 3.28 kPa, 2.74 kPa, and 3.56 kPa. The standard deviations of the estimated shear
elasticity values within the contours for these five measurements are 0.81 kPa, 0.84 kPa,
1.77 kPa, 1.33 kPa, and 1.60 kPa. These are also much larger than the standard deviations
obtained with the new model-based approach. Thus, the estimated shear elasticity values
estimated with the TOF method consistently have a larger variance for each measurement.
110

Measurement
Measurement
Measurement
Measurement
Measurement

1
2
3
4
5

Push beam focused at 35 mm
7th rib-cage gap 8th rib-cage gap
0.90 ± 0.29
1.24 ± 0.52
0.95 ± 0.30
1.01 ± 0.37
0.86 ± 0.29
1.00 ± 0.21
0.93 ± 0.31
1.07 ± 0.43
0.91 ± 0.32
1.12 ± 0.40

Push beam focused at 50 mm
7th rib-cage gap 8th rib-cage gap
1.36 ± 0.38
0.95 ± 0.36
1.06 ± 0.34
0.94 ± 0.28
1.25 ± 0.46
1.26 ± 0.49
1.41 ± 0.42
1.07 ± 0.42
1.47 ± 0.62
1.15 ± 0.38

Table 5.3: The average and the standard deviations of the shear viscosity values estimated
with the new model-based approach for 20 different measurements. The units are Pa·s for
each entry in the table.
Furthermore, for different measurement locations and for different push beams, significant
differences are observed in the estimated values.
Table 5.3 shows the estimated shear viscosity values provided by the new model-based
approach for the same twenty measurements. For the five shear wave measurements in the 7th
rib cage gap with a push beam focused at 35 mm, the average estimated shear viscosity values
are 0.90 Pa·s, 0.95 Pa·s, 0.86 Pa·s, 0.93 Pa·s, and 0.91 Pa·s. The standard deviations of the
estimated shear viscosity values within the contours for these five measurements are 0.29 Pa·s,
0.30 Pa·s, 0.29 Pa·s, 0.31 Pa·s, and 0.32 Pa·s. For the five shear wave measurements in the 8th
rib cage gap with a push beam focused at 35 mm, the average estimated shear viscosity values
are 1.24 Pa·s, 1.01 Pa·s, 1.00 Pa·s, 1.07 Pa·s, and 1.12 Pa·s. The standard deviations of the
estimated shear viscosity values within the contours for these five measurements are 0.52 Pa·s,
0.37 Pa·s, 0.21 Pa·s, 0.43 Pa·s, and 0.40 Pa·s. For the five shear wave measurements in the 7th
rib cage gap with a push beam focused at 50 mm, the average estimated shear viscosity values
are 1.36 Pa·s, 1.06 Pa·s, 1.25 Pa·s, 1.41 Pa·s, and 1.47 Pa·s. The standard deviations of the
estimated shear viscosity values within the contours for these five measurements are 0.38 Pa·s,
0.34 Pa·s, 0.46 Pa·s, 0.42 Pa·s, and 0.62 Pa·s. For the five shear wave measurements in the
8th rib cage gap with a push beam focused at 50 mm, the average estimated shear viscosity
values are 0.95 Pa·s, 0.94 Pa·s, 1.26 Pa·s, 1.07 Pa·s, and 1.15 Pa·s. The standard deviations
of the estimated shear viscosity values within the contours for these five measurements are
0.36 Pa·s, 0.28 Pa·s, 0.49 Pa·s, 0.42 Pa·s, and 0.38 Pa·s. The above results indicate that the
111

shear viscosity values estimated with the new model-based approach have small variances.
However, for different measurement locations and for different push beams, some differences
are evident in the estimated shear viscosity values.

5.5
5.5.1

Discussion
Shear waves measured with different push beams

Fig. 5.5 contains a comparison between measured shear wave particle velocities induced by
different push beams. In Figs. 5.5(a)-5.5(c), the maximum value of the shear wave particle
velocity in the x direction appears at x = 3.35 mm, x = 5.72 mm, and x = 8.08 mm at 2 ms,
4 ms, and 6 ms, respectively. This indicates that the shear waves consistently propagate away
from the push beam in the lateral direction as time progresses. In addition, the maximum
value in the z direction appears at z = 31.54 mm, z = 30.36 mm, and z = 28.39 mm at 2 ms,
4 ms, and 6 ms, respectively. This indicates that shear wave peak moves slightly upward as
time progresses. Also, the depths of the peak values of the shear wave particle velocities at all
three different time points are shallower than the depth of the focal peak of the push beam,
which is 35 mm. Similar behaviors are also observed in Figs. 5.5(d)-5.5(f). The maximum
values of the shear wave particle velocity in the x direction occur at x = 3.35 mm, x = 5.32
mm, and x = 8.48 mm at 2 ms, 4 ms, and 6 ms, respectively, and at z = 38.24 mm, z = 38.64
mm, and z = 35.09 mm at 2 ms, 4 ms, and 6 ms in the z direction, respectively. The depths
of the peak values for the shear wave particle velocity at all three times are also shallower
than the depth of the focal peak of the push beam, which is 50 mm. These effects occur due
to the diffraction of the shear wave. Phase aberrations and sound speed differences also play
a role, as different tissue layers and changes in sound speed also tend to shift the focal peak
of the push beam. Figs. 5.5(a)-5.5(c) and Figs. 5.5(d)-5.5(f) demonstrate that the FWHM
of the shear waves near the focal peak is different for shear waves induced by different push
beams. In Figs. 5.5(a)-5.5(c), the FWHM values of the shear wave particle velocity along the
z direction at the maximum position in x are 14.09 mm, 15.45 mm, and 15.21 mm for the

112

push beam focused at 35 mm, and the FWHM of the shear wave particle velocity induced
along the z direction at the maximum position in x are 24.24 mm, 24.12 mm, and 22.12 mm
for the push beam focused at 50 mm. This indicates that the FWHM values of the shear
wave particle velocities induced by the push beam focused at 50 mm are significantly larger
than the FWHM of the shear waves induced by the push beam focused at 35 mm. This is
due to the difference in the shape of the push beams as shown in Fig. 5.4, since the axial
extent of the push beam focused at 50 mm is significantly longer than the axial extent of
the push beam focused at 35 mm.

5.5.2

Shear wave parameter estimates in different locations

Fig. 5.6 contains the comparison between the measured and simulated shear wave particle
velocities with µ and ηs values estimated with the new model-based approach at different
locations in space. The results in Figs. 5.6(a) and 5.6(b) show that, for the shear wave
particle velocities measured at (x, z) = (5 mm, 30 mm) and at (x, z) = (6 mm, 35 mm),
the percent differences between the measured and simulated shear wave are 29.50% and
44.11%, respectively. Furthermore, in both of these locations, the measured and simulated
shear waves demonstrate good visual alignment. Also, in these two locations, the percent
difference is relatively small. For the shear wave particle velocities in Fig. 5.6(c) measured
at (x, z) = (3 mm, 30 mm) and in Fig. 5.6(d) measured and at (x, z) = (15 mm, 30 mm),
the percent difference between the measured and simulated shear wave particle velocities are
72.30% and 138.11%, respectively. Both of these results demonstrate poor visual alignment.
Thus, if the shear waves are measured too close or too far away from the push beam, poor
alignment between the measured and simulated shear wave particle velocities is expected
with large values for the percent difference.

113

5.5.3

Shear wave parameter estimates for in vivo liver

Tables 5.1 - 5.3 show the shear elasticity and shear viscosity values estimated with the TOF
method and the new model-based approach. The results shown in Table 5.1 reveal that
the standard deviations of the shear elasticity values estimated with the new model-based
approach within the contours for all test cases are quite small, as the smallest standard
deviation is 0.20 kPa, and the largest standard deviation is 0.49 kPa. This indicates that the
estimated shear elasticity values estimated with the new model-based approach have a small
variance within the contours for each measurement. However, for different measurement
locations, the shear elasticity estimates demonstrate some differences that are probably due
to inhomogeneities. The results shown in Table 5.2 reveal that the standard deviations of
the estimated shear elasticity values estimated in the contours with the TOF method are
sometimes very high, as the smallest standard deviation is 0.74 kPa, and the largest standard deviation is 6.89 kPa. These values are much larger than the standard deviations for
the shear elasticity values estimated with the new model-based approach. This indicates
that the shear elasticity values estimated with the TOF method have a larger variance for
each measurement compared with the new model-based approach. In addition, for different
measurement locations and for different push beams, some differences in the shear elasticity
estimates are evident. Table 5.3 reveals that the standard deviations of the estimated shear
viscosity values estimated in the contours with the new model-based approach are also small,
as the smallest standard deviation is 0.21 Pa·s, and the largest standard deviation is 0.62
Pa·s. This also suggests that the shear viscosity values estimated with the new model-based
approach have a small variance for each measurement within the contours. Overall, the
results indicate that the new model-based approach provides shear elasticity and shear viscosity estimates with smaller standard deviations and smaller percent differences compared
with the TOF method.

114

5.6

Conclusion

A new model-based approach that estimates the shear elasticity and shear viscosity with
measured shear wave particle velocities for in vivo human liver is introduced. In this new
approach, shear wave particle velocities are first calculated using an approximate viscoelastic
Green’s function model on a GPU, and the shear elasticity and shear viscosity values are then
estimated using a model-based approach that minimizes the difference between measured
and simulated shear wave particle velocities. Shear elasticity and shear viscosity images are
generated with this approach for all of the spatial points in a 2D FOV. The new approach is
applied to shear wave particle velocities obtained from in vivo human livers, and the results
show that this new approach successfully generates shear elasticity and shear viscosity images
for the in vivo human liver data. The results also indicate that the shear elasticity values
estimated with this approach are significantly smaller than the values estimated with the
TOF method, and the new approach yields more consistent results than the TOF method.

115

CHAPTER 6
CONCLUSION AND FUTURE WORK

6.1

Summary of the thesis

In this thesis, a new model-based approach that provides shear elasticity and shear viscosity estimates using shear wave particle velocities measured in viscoelastic media is introduced. In the new approach, the shear wave particle velocity is first calculated using a
viscoelastic model on a GPU. Then, the parameter combination that minimizes the percent
difference between the measured and simulated shear wave particle velocities is determined.
Shear wave parameters are estimated on a point-by-point basis to generate shear elasticity
and shear viscosity images. The new approach is applied to numerical simulations, viscoelastic shear wave phantoms, and in vivo human liver. The results show that this new
model-based approach successfully generates shear elasticity and shear viscosity images in
different viscoelastic media.
In chapter 3, shear waves induced by an acoustic radiation force are evaluated using
Green’s functions for simulations with elastic and viscoelastic models. The 3D acoustic radiation force is first simulated using the fast nearfield method and the angular spectrum
approach in FOCUS, and then shear waves are calculated by superposing the contributions from all of the source points in the 3D push beam using Green’s functions. Next,
these simulations of shear waves in elastic and viscoelastic models are accelerated using a
high-performance GPU. Diffraction effects are observed in the simulated shear wave particle
displacements in both elastic and viscoelastic models, while additional dispersion and attenuation effects are produced by the shear viscosity. Comparisons between the computation
times for serial and parallel implementations show that the nVidia Tesla K20 GPU accelerates these Green’s function calculations by factors of 45 and 700 for shear wave simulations
in elastic and viscoelastic models, respectively, relative to serial calculations on a desktop

116

computer with an Intel(R) Core(TM) i7-2600 3.40GHz CPU. When these Green’s functions
are evaluated for push beams with different spatial samplings, the results show that the
percent difference is larger in shear wave particle displacements calculated for push beams
with lower f/#. The results indicate that the percent difference usually decreases monotonically as the spatial sampling decreases, with the exception of the simulated shear waves in
a viscoelastic medium using f/1.5 and f/2 push beams, where the percent difference for both
of these reach a local minimum when the spatial sampling is λ/2. The results also indicate
that effective results are obtained using push beams sampled at λ/2, which yield percent differences less than 5%. Measured and simulated shear waves in elastic and viscoelastic shear
wave phantoms are compared, and the results suggest that the Green’s function simulation
approach provides an effective model for shear waves in elastic and viscoelastic shear wave
phantoms.
In chapter 4, a new model-based approach is introduced that provides estimates for the
shear elasticity and shear viscosity in viscoelastic shear wave phantoms. The shear wave
particle velocity is first calculated with different combinations of the shear elasticity and
shear viscosity on a GPU, then a new model-based approach is employed to determine the
parameter combination that yields the smallest percent difference between the measured and
simulated shear waves. The shear wave parameter estimates provided by the new modelbased approach are compared with the TOF method in simulation models, and the results
show that the new model-based approach provide more accurate result than the TOF method
when applied to simulation data obtained from a 3D viscoelastic shear wave model. Shear
wave parameters are then estimated in measured shear wave data collected from three different phantoms for three different push beams. The results suggest that the shear elasticity
values estimated with the TOF method are significantly higher than the shear elasticity values estimated with the new model-based approach, and the standard deviation is also higher
with the TOF approach. The results also indicate that the shear elasticity values estimated
with the new model-based approach consistently demonstrate much smaller standard devi-

117

ations for different push beams relative to those produced by the TOF method.
In chapter 5, the new model-based approach is accelerated with an approximate viscoelastic Green’s function model, and the approach is evaluated on shear wave particle velocities
obtained from in vivo human livers. Instead of calculating shear waves for combinations
of different shear elasticities and shear viscosities, shear waves are calculated with different
shear elasticities on GPU, and then the result is convolved with a loss model. This achieves
a dramatic reduction in the computation time. The shear elasticity and shear viscosity values are then estimated using a model-based approach that minimizes the difference between
measured and simulated shear wave particle velocities. Shear elasticity and shear viscosity
images are generated by evaluating this approach for all of the spatial points in a 2D FOV.
The new approach is applied to shear wave particle velocities obtained from in vivo human
livers, and the results show that this new approach successfully generates shear elasticity
and shear viscosity images for the in vivo human livers. The results also indicate that the
shear elasticity values estimated with this approach are significantly smaller than the values
estimated with the TOF method, and the new approach yields more consistent results than
the TOF method.

6.2

Future work

One possible future direction is to develop a more advanced model to represent viscous
loss. In chapter 3 - 5, Voigt’s model describes the viscoelasticity of soft tissue. The Voigt
model, which is widely used, provides an effective description of some viscoelastic media.
However, this model does not describe shear waves in soft tissue exactly. Measurements
suggest that the attenuation of shear waves in soft tissue follows a power law with an exponent
between 0 and 2. Power law attenuation of shear waves is observed in the brain [181, 182],
breast [183], muscle [182], liver [181,182,184,185], prostate [186,187], and arterial tissue [188].
Thus, fractional calculus models have great potential for modeling shear waves propagating
in soft tissue.

118

Another possible direction for future work is to investigate different models for diffracting shear waves. In chapter 5, the diffraction of shear waves is calculated by evaluating the
shear waves induced by a 3D acoustic radiation force in elastic media. Although this is faster
than evaluating shear waves with combinations of all µ and ηs values, this calculation is still
quite time consuming. As an alternative, line source models could be used to describe the
diffraction of shear waves near the focal peak of the push beam, since shear wave propagation at the depth of the focal peak of the acoustic radiation force is similar to that of a
cylindrical wave. The line source model can potentially provide an effective approximation
for the diffraction of shear waves, where the computation time could be reduced if such an
approximation could be made.
In addition, more work could be done to implement the estimation procedure on a GPU.
Currently, only the shear wave simulations are performed in parallel on a GPU. If the new
model-based approach is performed on a GPU, this should increase the speed of the estimation procedure significantly. Also, the new estimation approach should be applied to shear
waves measured in other organs, such as the kidney, thyroid, and prostate.

119

BIBLIOGRAPHY

120

BIBLIOGRAPHY

[1]

F. Adams et al., The Genuine Works of Hippocrates, vol. 1. Sydenham society, 1849.

[2]

J. N. Wolfe, “Xeroradiography of the breast,” Oncology, vol. 23, no. 2, pp. 113–119,
1969.

[3]

S. Weiss, “Liver deaths and their prevention: How the danger can be recognized and
avoided by use of preoperative and postoperative diagnostic measures,” The American
Journal of Surgery, vol. 23, no. 1, pp. 96–101, 1934.

[4]

J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imaging elastic
properties of biological tissues,” Annual Review of Biomedical Engineering, vol. 5,
no. 1, pp. 57–78, 2003.

[5]

K. J. Parker, L. S. Taylor, S. Gracewski, and D. J. Rubens, “A unified view of imaging
the elastic properties of tissue,” The Journal of the Acoustical Society of America,
vol. 117, no. 5, pp. 2705–2712, 2005.

[6]

J. Ophir, S. Srinivasan, R. Righetti, and A. Thittai, “Elastography: A decade of
progress (2000-2010),” Current Medical Imaging Reviews, vol. 7, no. 4, pp. 292–312,
2011.

[7]

K. Parker, M. Doyley, and D. Rubens, “Imaging the elastic properties of tissue: the
20 year perspective,” Physics in Medicine and Biology, vol. 56, no. 1, p. R1, 2011.

[8]

M. L. Palmeri and K. R. Nightingale, “Acoustic radiation force-based elasticity imaging
methods,” Interface Focus, vol. 1, no. 4, pp. 553–564, 2011.

[9]

A. Sarvazyan, T. J. Hall, M. W. Urban, M. Fatemi, S. R. Aglyamov, and B. S. Garra,
“An overview of elastography–an emerging branch of medical imaging,” Current Medical Imaging Reviews, vol. 7, no. 4, p. 255, 2011.

[10] J. R. Doherty, G. E. Trahey, K. R. Nightingale, and M. L. Palmeri, “Acoustic radiation
force elasticity imaging in diagnostic ultrasound,” IEEE Transactions on Ultrasonics,
Ferroelectrics and Frequency Control, vol. 60, no. 4, pp. 685–701, 2013.
[11] C. Kasai, K. Namekawa, A. Koyano, and R. Omoto, “Real-time two-dimensional blood
flow imaging using an autocorrelation technique,” IEEE Transactions on Sonics and
Ultrasonics, vol. 32, no. 3, pp. 458–464, 1985.
[12] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: a quantitative method for imaging the elasticity of biological tissues,” Ultrasonic Imaging,
vol. 13, no. 2, pp. 111–134, 1991.
[13] I. Cespedes, J. Ophir, H. Ponnekanti, and N. Maklad, “Elastography: elasticity imaging using ultrasound with application to muscle and breast in vivo,” Ultrasonic Imaging, vol. 15, no. 2, pp. 73–88, 1993.
121

[14] M. Niitsu, A. Michizaki, A. Endo, H. Takei, and O. Yanagisawa, “Muscle hardness
measurement by using ultrasound elastography: a feasibility study,” Acta Radiologica,
vol. 52, no. 1, pp. 99–105, 2011.
[15] O. Yanagisawa, M. Niitsu, T. Kurihara, and T. Fukubayashi, “Evaluation of human
muscle hardness after dynamic exercise with ultrasound real-time tissue elastography:
a feasibility study,” Clinical Radiology, vol. 66, no. 9, pp. 815–819, 2011.
[16] S. Chan, P. Fung, N. Ng, T. Ngan, M. Chong, C. Tang, J. He, and Y. Zheng, “Dynamic
changes of elasticity, cross-sectional area, and fat infiltration of multifidus at different
postures in men with chronic low back pain,” The Spine Journal, vol. 12, no. 5, pp. 381–
388, 2012.
[17] B. S. Garra, E. I. Cespedes, J. Ophir, S. R. Spratt, R. A. Zuurbier, C. M. Magnant, and
M. F. Pennanen, “Elastography of breast lesions: initial clinical results.,” Radiology,
vol. 202, no. 1, pp. 79–86, 1997.
[18] A. Itoh, E. Ueno, E. Tohno, H. Kamma, H. Takahashi, T. Shiina, M. Yamakawa, and
T. Matsumura, “Breast disease: clinical application of US elastography for diagnosis,”
Radiology, vol. 239, no. 2, pp. 341–350, 2006.
[19] A. Thitaikumar, L. M. Mobbs, C. M. Kraemer-Chant, B. S. Garra, and J. Ophir,
“Breast tumor classification using axial shear strain elastography: a feasibility study,”
Physics in Medicine and Biology, vol. 53, no. 17, p. 4809, 2008.
[20] W. Moon, S. Chang, C. Huang, and R. Chang, “Breast tumor classification using fuzzy
clustering for breast elastography,” Ultrasound in Medicine & Biology, vol. 37, no. 5,
pp. 700–708, 2011.
[21] A. Lyshchik, T. Higashi, R. Asato, S. Tanaka, J. Ito, J. J. Mai, C. Pellot-Barakat,
M. F. Insana, A. B. Brill, T. Saga, M. Hiraoka, and T. Togashi, “Thyroid gland tumor
diagnosis at US elastography,” Radiology, vol. 237, no. 1, pp. 202–211, 2005.
[22] T. Rago, F. Santini, M. Scutari, A. Pinchera, and P. Vitti, “Elastography: new developments in ultrasound for predicting malignancy in thyroid nodules,” The Journal of
Clinical Endocrinology & Metabolism, vol. 92, no. 8, pp. 2917–2922, 2007.
[23] C. Asteria, A. Giovanardi, A. Pizzocaro, L. Cozzaglio, A. Morabito, F. Somalvico,
and A. Zoppo, “US-elastography in the differential diagnosis of benign and malignant
thyroid nodules,” Thyroid, vol. 18, no. 5, pp. 523–531, 2008.
[24] R. H. Cochlin, D. L.and Ganatra and D. F. R. Griffiths, “Elastography in the detection
of prostatic cancer,” Clinical Radiology, vol. 57, no. 11, pp. 1014–1020, 2002.
[25] K. König, U. Scheipers, A. Pesavento, A. Lorenz, H. Ermert, and T. Senge, “Initial
experiences with real-time elastography guided biopsies of the prostate,” The Journal
of Urology, vol. 174, no. 1, pp. 115–117, 2005.

122

[26] G. Salomon, J. Köllerman, I. Thederan, F. K. H. Chun, L. Budäus, T. Schlomm, H. Isbarn, H. Heinzer, H. Huland, and M. Graefen, “Evaluation of prostate cancer detection
with ultrasound real-time elastography: a comparison with step section pathological
analysis after radical prostatectomy,” European Urology, vol. 54, no. 6, pp. 1354–1362,
2008.
[27] R. Righetti, F. Kallel, R. J. Stafford, R. E. Price, T. A. Krouskop, J. D. Hazle, and
J. Ophir, “Elastographic characterization of HIFU-induced lesions in canine livers,”
Ultrasound in Medicine & Biology, vol. 25, no. 7, pp. 1099–1113, 1999.
[28] T. Varghese, J. A. Zagzebski, and F. T. Lee, “Elastographic imaging of thermal lesions
in the liver in vivo following radiofrequency ablation: preliminary results,” Ultrasound
in Medicine & Biology, vol. 28, no. 11, pp. 1467–1473, 2002.
[29] T. Varghese, U. Techavipoo, W. Liu, J. A. Zagzebski, Q. Chen, G. Frank, and F. T.
Lee Jr, “Elastographic measurement of the area and volume of thermal lesions resulting from radiofrequency ablation: pathologic correlation,” American Journal of
Roentgenology, vol. 181, no. 3, pp. 701–707, 2003.
[30] C. L. de Korte, A. F. W. van der Steen, E. I. Céspedes, G. Pasterkamp, S. G. Carlier,
F. Mastik, A. H. Schoneveld, P. W. Serruys, and N. Bom, “Characterization of plaque
components and vulnerability with intravascular ultrasound elastography,” Physics in
Medicine and Biology, vol. 45, no. 6, p. 1465, 2000.
[31] A. P. Sarvazyan, S. Emelianov, and A. R. Skovoroda, “Intracavity ultrasonic device
for elasticity imaging,” Nov. 30 1993. US Patent 5,265,612.
[32] C. L. de Korte, M. J. Sierevogel, F. Mastik, C. Strijder, J. A. Schaar, E. Velema,
G. Pasterkamp, P. W. Serruys, and A. F. W. van der Steen, “Identification of
atherosclerotic plaque components with intravascular ultrasound elastography in vivo,”
Circulation, vol. 105, no. 14, pp. 1627–1630, 2002.
[33] J. J. Mai and M. F. Insana, “Strain imaging of internal deformation,” Ultrasound in
Medicine & Biology, vol. 28, no. 11, pp. 1475–1484, 2002.
[34] R. L. Maurice, J. Ohayon, Y. Frétigny, M. Bertrand, G. Soulez, and G. Cloutier,
“Noninvasive vascular elastography: Theoretical framework,” IEEE Transactions on
Medical Imaging, vol. 23, no. 2, pp. 164–180, 2004.
[35] S. Korukonda, R. Nayak, N. Carson, G. Schifitto, V. Dogra, and M. M. Doyley, “Noninvasive vascular elastography using plane-wave and sparse-array imaging,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 60, no. 2, pp. 332–
342, 2013.
[36] K. G. Sabra, S. Conti, P. Roux, and W. A. Kuperman, “Passive in vivo elastography
from skeletal muscle noise,” Applied Physics Letters, vol. 90, no. 19, p. 194101, 2007.

123

[37] W. Kuo, D. Jian, T. Wang, and Y.-C. Wang, “Neck muscle stiffness quantified by
sonoelastography is correlated with body mass index and chronic neck pain symptoms,”
Ultrasound in Medicine & Biology, vol. 39, no. 8, pp. 1356–1361, 2013.
[38] E. E. Konofagou, J. Dhooge, and J. Ophir, “Myocardial elastographya feasibility study
in vivo,” Ultrasound in Medicine & Biology, vol. 28, no. 4, pp. 475–482, 2002.
[39] T. Varghese, J. A. Zagzebski, P. Rahko, and C. S. Breburda, “Ultrasonic imaging
of myocardial strain using cardiac elastography,” Ultrasonic Imaging, vol. 25, no. 1,
pp. 1–16, 2003.
[40] H. Kanai, “Propagation of spontaneously actuated pulsive vibration in human heart
wall and in vivo viscoelasticity estimation,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 52, no. 11, pp. 1931–1942, 2005.
[41] K. R. Nightingale, M. L. Palmeri, R. W. Nightingale, and G. E. Trahey, “On the feasibility of remote palpation using acoustic radiation force,” The Journal of the Acoustical
Society of America, vol. 110, no. 1, pp. 625–634, 2001.
[42] A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov,
“Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics,”
Ultrasound in Medicine & Biology, vol. 24, no. 9, pp. 1419–1435, 1998.
[43] K. Nightingale, M. S. Soo, R. Nightingale, and G. Trahey, “Acoustic radiation force
impulse imaging: in vivo demonstration of clinical feasibility,” Ultrasound in Medicine
& Biology, vol. 28, no. 2, pp. 227–235, 2002.
[44] C. Fierbinteanu-Braticevici, D. Andronescu, R. Usvat, D. Cretoiu, C. Baicus, and
G. Marinoschi, “Acoustic radiation force imaging sonoelastography for noninvasive
staging of liver fibrosis,” World Journal of Gastroenterology, vol. 15, no. 44, p. 5525,
2009.
[45] B. Fahey, R. Nelson, D. Bradway, S. Hsu, D. Dumont, and G. Trahey, “In vivo visualization of abdominal malignancies with acoustic radiation force elastography,” Physics
in Medicine and Biology, vol. 53, no. 1, p. 279, 2008.
[46] K. S. S. Bhatia, D. P. Rasalkar, Y. P. Lee, K. T. Wong, A. D. King, H. Y. Yuen,
and A. T. Ahuja, “Cystic change in thyroid nodules: a confounding factor for realtime qualitative thyroid ultrasound elastography,” Clinical Radiology, vol. 66, no. 9,
pp. 799–807, 2011.
[47] Y. Zhang, Y. He, H. Xu, X. Xu, C. Liu, L. Guo, L. Liu, and J. Xu, “Virtual touch
tissue imaging on acoustic radiation force impulse elastography,” Journal of Ultrasound
in Medicine, vol. 33, no. 4, pp. 585–595, 2014.
[48] A. C. Sharma, M. S. Soo, G. E. Trahey, and K. R. Nightingale, “Acoustic radiation
force impulse imaging of in vivo breast masses,” in 2004 IEEE Ultrasonics Symposium,
vol. 1, pp. 728–731, IEEE, 2004.

124

[49] M. Tozaki, S. Isobe, M. Yamaguchi, Y. Ogawa, K. Homma, M. Saito, C. Joo, and
E. Fukuma, “Ultrasonographic elastography of the breast using acoustic radiation force
impulse technology: preliminary study,” Japanese Journal of Radiology, vol. 29, no. 6,
pp. 452–456, 2011.
[50] L. Zhai, J. Madden, W.-C. Foo, M. L. Palmeri, V. Mouraviev, T. J. Polascik, and
K. R. Nightingale, “Acoustic radiation force impulse imaging of human prostates ex
vivo,” Ultrasound in Medicine & Biology, vol. 36, no. 4, pp. 576–588, 2010.
[51] L. Zhai, T. J. Polascik, W.-C. Foo, S. Rosenzweig, M. L. Palmeri, J. Madden, and K. R.
Nightingale, “Acoustic radiation force impulse imaging of human prostates: Initial in
vivo demonstration,” Ultrasound in Medicine & Biology, vol. 38, no. 1, pp. 50–61,
2012.
[52] J. Allen, D. Dumont, B. Fahey, E. Miller, J. Dahl, and G. Trahey, “Lower-limb vascular imaging with acoustic radiation force elastography: demonstration of in vivo
feasibility,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control,
vol. 56, no. 5, pp. 931–944, 2009.
[53] D. Dumont, R. H. Behler, T. C. Nichols, E. P. Merricks, and C. M. Gallippi, “ARFI
imaging for noninvasive material characterization of atherosclerosis,” Ultrasound in
Medicine & Biology, vol. 32, no. 11, pp. 1703–1711, 2006.
[54] G. E. Trahey, M. L. Palmeri, R. C. Bentley, and K. R. Nightingale, “Acoustic radiation
force impulse imaging of the mechanical properties of arteries: In vivo and ex vivo
results,” Ultrasound in Medicine & Biology, vol. 30, no. 9, pp. 1163–1171, 2004.
[55] J. J. Dahl, D. M. Dumont, J. D. Allen, E. M. Miller, and G. E. Trahey, “Acoustic radiation force impulse imaging for noninvasive characterization of carotid artery
atherosclerotic plaques: a feasibility study,” Ultrasound in Medicine & Biology, vol. 35,
no. 5, pp. 707–716, 2009.
[56] J. R. Doherty, D. M. Dumont, G. E. Trahey, and M. L. Palmeri, “Acoustic radiation force impulse imaging of vulnerable plaques: a finite element method parametric
analysis,” Journal of Biomechanics, vol. 46, no. 1, pp. 83–90, 2013.
[57] S. J. Hsu, R. R. Bouchard, D. M. Dumont, P. D. Wolf, and G. E. Trahey, “In vivo
assessment of myocardial stiffness with acoustic radiation force impulse imaging,” Ultrasound in Medicine & Biology, vol. 33, no. 11, pp. 1706–1719, 2007.
[58] S. J. Hsu, J. L. Hubert, S. W. Smith, and G. E. Trahey, “Intracardiac echocardiography
and acoustic radiation force impulse imaging of a dynamic ex vivo ovine heart model,”
Ultrasonic Imaging, vol. 30, no. 2, pp. 63–77, 2008.
[59] S. A. Eyerly, S. J. Hsu, S. H. Agashe, G. E. Trahey, Y. Li, and P. D. Wolf, “An in
vitro assessment of acoustic radiation force impulse imaging for visualizing cardiac
radiofrequency ablation lesions,” Journal of Cardiovascular Electrophysiology, vol. 21,
no. 5, pp. 557–563, 2010.

125

[60] K. Nightingale, S. McAleavey, and G. Trahey, “Shear-wave generation using acoustic
radiation force: in vivo and ex vivo results,” Ultrasound in Medicine & Biology, vol. 29,
no. 12, pp. 1715–1723, 2003.
[61] M. L. Palmeri, M. H. Wang, J. J. Dahl, K. D. Frinkley, and K. R. Nightingale, “Quantifying hepatic shear modulus in vivo using acoustic radiation force,” Ultrasound in
Medicine & Biology, vol. 34, no. 4, pp. 546–558, 2008.
[62] M. L. Palmeri, M. H. Wang, N. C. Rouze, M. F. Abdelmalek, C. D. Guy, B. Moser,
A. M. Diehl, and K. R. Nightingale, “Noninvasive evaluation of hepatic fibrosis using
acoustic radiation force-based shear stiffness in patients with nonalcoholic fatty liver
disease,” Journal of Hepatology, vol. 55, no. 3, pp. 666–672, 2011.
[63] M. Friedrich-Rust, K. Wunder, S. Kriener, F. Sotoudeh, S. Richter, J. Bojunga, E. Herrmann, T. Poynard, C. F. Dietrich, J. Vermehren, et al., “Liver fibrosis in viral hepatitis: Noninvasive assessment with acoustic radiation force impulse imaging versus
transient elastography 1,” Radiology, vol. 252, no. 2, pp. 595–604, 2009.
[64] M. Haque, C. Robinson, D. Owen, E. M. Yoshida, and A. Harris, “Comparison of
acoustic radiation force impulse imaging (ARFI) to liver biopsy histologic scores in
the evaluation of chronic liver disease: A pilot study,” Ann Hepatol, vol. 9, no. 3,
pp. 289–93, 2010.
[65] S. Horster, P. Mandel, R. Zachoval, and D. Clevert, “Comparing acoustic radiation
force impulse imaging to transient elastography to assess liver stiffness in healthy volunteers with and without valsalva manoeuvre,” Clinical Hemorheology and Microcirculation, vol. 46, no. 2, pp. 159–168, 2010.
[66] M. Lupsor, R. Badea, H. Stefanescu, Z. Sparchez, H. Branda, A. Serban, and A. Maniu, “Performance of a new elastographic method (ARFI technology) compared to
unidimensional transient elastography in the noninvasive assessment of chronic hepatitis c. preliminary results.,” Journal of Gastrointestinal & Liver Diseases, vol. 18,
no. 3, 2009.
[67] F. Sebag, J. Vaillant-Lombard, J. Berbis, V. Griset, J. F. Henry, P. Petit, and C. Oliver,
“Shear wave elastography: a new ultrasound imaging mode for the differential diagnosis
of benign and malignant thyroid nodules,” The Journal of Clinical Endocrinology &
Metabolism, vol. 95, no. 12, pp. 5281–5288, 2010.
[68] R. Goertz, K. Amann, R. Heide, T. Bernatik, M. Neurath, and D. Strobel, “An abdominal and thyroid status with acoustic radiation force impulse elastometry–a feasibility
study: Acoustic radiation force impulse elastometry of human organs,” European Journal of Radiology, vol. 80, no. 3, pp. e226–e230, 2011.
[69] V. Ianculescu, L. M. Ciolovan, A. Dunant, P. Vielh, C. Mazouni, S. Delaloge, C. Dromain, A. Blidaru, and C. Balleyguier, “Added value of virtual touch IQ shear wave
elastography in the ultrasound assessment of breast lesions,” European Journal of Radiology, vol. 83, no. 5, pp. 773–777, 2014.
126

[70] C. Balleyguier, S. Canale, W. B. Hassen, P. Vielh, E. H. Bayou, M. C. Mathieu,
C. Uzan, C. Bourgier, and C. Dromain, “Breast elasticity: principles, technique, results: an update and overview of commercially available software,” European Journal
of Radiology, vol. 82, no. 3, pp. 427–434, 2013.
[71] K. F. Stock, B. S. Klein, M. T. Cong, C. Regenbogen, S. Kemmner, M. Büttner,
S. Wagenpfeil, E. Matevossian, L. Renders, U. Heemann, and C. Küchle, “ARFI-based
tissue elasticity quantification and kidney graft dysfunction: first clinical experiences,”
Clinical Hemorheology and Microcirculation, vol. 49, no. 1-4, pp. 527–535, 2011.
[72] X. Zheng, P. Ji, H. Mao, and J. Hu, “A comparison of virtual touch tissue quantification
and digital rectal examination for discrimination between prostate cancer and benign
prostatic hyperplasia,” Radiology and Oncology, vol. 46, no. 1, pp. 69–74, 2012.
[73] K. Boehm, G. Salomon, B. Beyer, J. Schiffmann, K. Simonis, M. Graefen, and L. Budaeus, “Shear wave elastography for localization of prostate cancer lesions and assessment of elasticity thresholds: implications for targeted biopsies and active surveillance
protocols,” The Journal of Urology, vol. 193, no. 3, pp. 794–800, 2015.
[74] 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.
[75] M. Tanter, J. Bercoff, L. Sandrin, and M. Fink, “Ultrafast compound imaging for 2-d
motion vector estimation: Application to transient elastography,” IEEE Transactions
on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 49, no. 10, pp. 1363–1374,
2002.
[76] M. Tanter and M. Fink, “Ultrafast imaging in biomedical ultrasound,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 61, no. 1, pp. 102–119,
2014.
[77] M. Muller, J.-L. Gennisson, T. Deffieux, M. Tanter, and M. Fink, “Quantitative viscoelasticity mapping of human liver using supersonic shear imaging: Preliminary in
vivo feasability study,” Ultrasound in Medicine & Biology, vol. 35, no. 2, pp. 219–229,
2009.
[78] M. Tanter, D. Touboul, J. Gennisson, J. Bercoff, and M. Fink, “High-resolution quantitative imaging of cornea elasticity using supersonic shear imaging,” IEEE Transactions
on Medical Imaging, vol. 28, no. 12, pp. 1881–1893, 2009.
[79] É. Bavu, J. Gennisson, M. Couade, J. Bercoff, V. Mallet, M. Fink, A. Badel, A. ValletPichard, B. Nalpas, M. Tanter, et al., “Noninvasive in vivo liver fibrosis evaluation
using supersonic shear imaging: a clinical study on 113 hepatitis c virus patients,”
Ultrasound in Medicine & Biology, vol. 37, no. 9, pp. 1361–1373, 2011.
[80] G. Ferraioli, C. Tinelli, M. Zicchetti, E. Above, G. Poma, M. Di Gregorio, and C. Filice, “Reproducibility of real-time shear wave elastography in the evaluation of liver
elasticity,” European Journal of Radiology, vol. 81, no. 11, pp. 3102–3106, 2012.
127

[81] M. Tanter, J. Bercoff, A. Athanasiou, T. Deffieux, J.-L. Gennisson, G. Montaldo,
M. Muller, A. Tardivon, and M. Fink, “Quantitative assessment of breast lesion viscoelasticity: initial clinical results using supersonic shear imaging,” Ultrasound in
Medicine & Biology, vol. 34, no. 9, pp. 1373–1386, 2008.
[82] A. Athanasiou, A. Tardivon, M. Tanter, B. Sigal-Zafrani, J. Bercoff, T. Deffieux,
J. Gennisson, M. Fink, and S. Neuenschwander, “Breast lesions: quantitative elastography with supersonic shear imagingpreliminary results,” Radiology, vol. 256, no. 1,
pp. 297–303, 2010.
[83] T. Deffieux, J.-L. Gennisson, M. Tanter, and M. Fink, “Assessment of the mechanical
properties of the musculoskeletal system using 2-d and 3-d very high frame rate ultrasound,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control,
vol. 55, no. 10, pp. 2177–2190, 2008.
[84] J. Gennisson, T. Deffieux, E. Macé, G. Montaldo, M. Fink, and M. Tanter, “Viscoelastic and anisotropic mechanical properties of in vivo muscle tissue assessed by supersonic
shear imaging,” Ultrasound in Medicine & Biology, vol. 36, no. 5, pp. 789–801, 2010.
[85] A. Nordez and F. Hug, “Muscle shear elastic modulus measured using supersonic
shear imaging is highly related to muscle activity level,” Journal of Applied Physiology,
vol. 108, no. 5, pp. 1389–1394, 2010.
[86] S. A. McAleavey, M. Menon, and J. Orszulak, “Shear-modulus estimation by application of spatially-modulated impulsive acoustic radiation force,” Ultrasonic Imaging,
vol. 29, no. 2, pp. 87–104, 2007.
[87] S. McAleavey, E. Collins, J. Kelly, E. Elegbe, and M. Menon, “Validation of SMURF
estimation of shear modulus in hydrogels,” Ultrasonic Imaging, vol. 31, no. 2, pp. 131–
150, 2009.
[88] P. Song, H. Zhao, A. Manduca, M. W. Urban, J. F. Greenleaf, and S. Chen, “Combpush ultrasound shear elastography (CUSE): a novel method for two-dimensional shear
elasticity imaging of soft tissues,” IEEE Transactions on Medical Imaging, vol. 31,
no. 9, pp. 1821–1832, 2012.
[89] S. Chen, H. Zhao, P. Song, M. W. Urban, and J. Greenleaf, “Shear wave speed measurement using an unfocused ultrasound push beam,” in 2011 IEEE International
Ultrasonics Symposium, pp. 1171–1174, IEEE, 2011.
[90] P. Song, M. W. Urban, A. Manduca, H. Zhao, J. F. Greenleaf, and S. Chen, “Combpush ultrasound shear elastography (CUSE) with various ultrasound push beams,”
IEEE Transactions on Medical Imaging, vol. 32, no. 8, pp. 1435–1447, 2013.
[91] M. Denis, M. Bayat, M. Mehrmohammadi, A. Gregory, P. Song, D. H. Whaley,
S. Pruthi, S. Chen, M. Fatemi, and A. Alizad, “Update on breast cancer detection
using comb-push ultrasound shear elastography,” IEEE Transactions on Ultrasonics,
Ferroelectrics, and Frequency Control, vol. 62, no. 9, pp. 1644–1650, 2015.

128

[92] P. Song, D. C. Mellema, S. P. Sheedy, D. D. Meixner, R. M. Karshen, M. W. Urban,
A. Manduca, W. Sanchez, M. R. Callstrom, J. F. Greenleaf, and S. Chen, “Performance
of 2-dimensional ultrasound shear wave elastography in liver fibrosis detection using
magnetic resonance elastography as the reference standard,” Journal of Ultrasound in
Medicine, vol. 35, no. 2, pp. 401–412, 2016.
[93] R. M. Lerner, S. Huang, and K. J. Parker, “sonoelasticity images derived from ultrasound signals in mechanically vibrated tissues,” Ultrasound in Medicine & Biology,
vol. 16, no. 3, pp. 231–239, 1990.
[94] M. Fatemi and J. F. Greenleaf, “Ultrasound-stimulated vibro-acoustic spectrography,”
Science, vol. 280, no. 5360, pp. 82–85, 1998.
[95] M. W. Urban, M. Fatemi, and J. F. Greenleaf, “Modulation of ultrasound to produce
multifrequency radiation force),” The Journal of the Acoustical Society of America,
vol. 127, no. 3, pp. 1228–1238, 2010.
[96] K. J. Parker, S. R. Huang, R. A. Musulin, and R. M. Lerner, “Tissue response to
mechanical vibrations for sonoelasticity imaging,” Ultrasound in Medicine & Biology,
vol. 16, no. 3, pp. 241–246, 1990.
[97] S. F. Levinson, M. Shinagawa, and T. Sato, “Sonoelastic determination of human
skeletal muscle elasticity,” Journal of Biomechanics, vol. 28, no. 10, pp. 1145–1154,
1995.
[98] K. J. Parker and R. M. Lerner, “Sonoelasticity of organs: shear waves ring a bell,”
Journal of Ultrasound in Medicine, vol. 11, no. 8, pp. 387–392, 1992.
[99] M. Sanada, M. Ebara, H. Fukuda, M. Yoshikawa, N. Sugiura, H. Saisho, Y. Yamakoshi,
K. Ohmura, A. Kobayashi, and F. Kondoh, “Clinical evaluation of sonoelasticity measurement in liver using ultrasonic imaging of internal forced low-frequency vibration,”
Ultrasound in Medicine & Biology, vol. 26, no. 9, pp. 1455–1460, 2000.
[100] S. Y. Emelianov, M. A. Lubinski, W. F. Weitzel, R. C. Wiggins, A. R. Skovoroda, and
M. O’donnell, “Elasticity imaging for early detection of renal pathology,” Ultrasound
in Medicine & Biology, vol. 21, no. 7, pp. 871–883, 1995.
[101] K. J. Parker, R. M. Lerner, and S. Huang, “Method and apparatus for breast imaging and tumor detection using modal vibration analysis,” Mar. 31 1992. US Patent
5,099,848.
[102] D. J. Rubens, M. A. Hadley, S. K. Alam, L. Gao, R. D. Mayer, and K. J. Parker,
“Sonoelasticity imaging of prostate cancer: in vitro results.,” Radiology, vol. 195, no. 2,
pp. 379–383, 1995.
[103] L. S. Taylor, D. J. Rubens, B. C. Porter, Z. Wu, R. B. Baggs, P. A. di Sant’Agnese,
G. Nadasdy, D. Pasternack, E. M. Messing, P. Nigwekar, and K. J. Parker, “Prostate
cancer: three-dimensional sonoelastography for in vitro detection,” Radiology, vol. 237,
no. 3, pp. 981–985, 2005.
129

[104] M. Fatemi and J. F. Greenleaf, “Vibro-acoustography: An imaging modality based
on ultrasound-stimulated acoustic emission,” Proceedings of the National Academy of
Sciences, vol. 96, no. 12, pp. 6603–6608, 1999.
[105] A. Alizad, L. E. Wold, J. F. Greenleaf, and M. Fatemi, “Imaging mass lesions by vibroacoustography: modeling and experiments,” IEEE Transactions on Medical Imaging,
vol. 23, no. 9, pp. 1087–1093, 2004.
[106] M. Fatemi, L. E. Wold, A. Alizad, and J. F. Greenleaf, “Vibro-acoustic tissue mammography,” IEEE Transactions on Medical Imaging, vol. 21, no. 1, pp. 1–8, 2002.
[107] A. Alizad, M. Fatemi, L. E. Wold, and J. F. Greenleaf, “Performance of vibroacoustography in detecting microcalcifications in excised human breast tissue: A study
of 74 tissue samples,” IEEE Transactions on Medical Imaging, vol. 23, no. 3, pp. 307–
312, 2004.
[108] A. Allizad, F. G. Mitri, R. R. Kinnick, J. F. Greenleaf, and M. Fatemi, “Vibroacoustography of thyroid,” The Journal of the Acoustical Society of America, vol. 122,
no. 5, pp. 3026–3026, 2007.
[109] A. Alizad, F. G. Mitri, R. R. Kinnick, J. F. Greenleaf, and M. Fatemi, “Detection
of thyroid nodules by a novel acoustic imaging method,” in International Congress of
Ultrasound, vol. 35, 2009.
[110] F. G. Mitri, B. J. Davis, A. Alizad, J. F. Greenleaf, T. M. Wilson, L. A. Mynderse, and
M. Fatemi, “Prostate cryotherapy monitoring using vibroacoustography: preliminary
results of an ex vivo study and technical feasibility,” IEEE Transactions on Biomedical
Engineering, vol. 55, no. 11, pp. 2584–2592, 2008.
[111] F. G. Mitri, B. J. Davis, M. W. Urban, A. Alizad, J. F. Greenleaf, G. H. Lischer,
T. M. Wilson, and M. Fatemi, “Vibro-acoustography imaging of permanent prostate
brachytherapy seeds in an excised human prostate–preliminary results and technical
feasibility,” Ultrasonics, vol. 49, no. 3, pp. 389–394, 2009.
[112] C. Pislaru, B. Kantor, R. R. Kinnick, J. L. Anderson, M.-C. Aubry, M. W. Urban,
M. Fatemi, and J. F. Greenleaf, “In vivo vibroacoustography of large peripheral arteries,” Investigative Radiology, vol. 43, no. 4, p. 243, 2008.
[113] E. E. Konofagou and K. Hynynen, “Localized harmonic motion imaging: theory, simulations and experiments,” Ultrasound in Medicine & Biology, vol. 29, no. 10, pp. 1405–
1413, 2003.
[114] C. Maleke and E. E. Konofagou, “Harmonic motion imaging for focused ultrasound
(HMIFU): a fully integrated technique for sonication and monitoring of thermal ablation in tissues,” Physics in Medicine and Biology, vol. 53, no. 6, p. 1773, 2008.
[115] S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity and viscosity from
measurement of shear wave speed dispersion,” The Journal of the Acoustical Society
of America, vol. 115, no. 6, pp. 2781–2785, 2004.
130

[116] S. Chen, M. W. Urban, C. Pislaru, R. Kinnick, Y. Zheng, A. Yao, and J. F. Greenleaf,
“Shear wave dispersion ultrasound vibrometry (SDUV) for measuring tissue elasticity and viscosity,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency
Control, vol. 56, no. 1, pp. 55–62, 2009.
[117] C. Amador, M. W. Urban, S. Chen, and J. F. Greenleaf, “Loss tangent and complex modulus estimated by acoustic radiation force creep and shear wave dispersion,”
Physics in Medicine and Biology, vol. 57, no. 5, p. 1263, 2012.
[118] M. W. Urban and J. F. Greenleaf, “A kramers–kronig-based quality factor for shear
wave propagation in soft tissue,” Physics in Medicine and Biology, vol. 54, no. 19,
p. 5919, 2009.
[119] M. Bernal, I. Nenadic, M. W. Urban, and J. F. Greenleaf, “Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating
modes,” The Journal of the Acoustical Society of America, vol. 129, no. 3, pp. 1344–
1354, 2011.
[120] C. Amador, M. W. Urban, L. V. Warner, and J. F. Greenleaf, “In vitro renal cortex
elasticity and viscosity measurements with shearwave dispersion ultrasound vibrometry
(SDUV) on swine kidney,” in Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pp. 4428–4431, IEEE, 2009.
[121] C. Amador, M. W. Urban, S. Chen, and J. F. Greenleaf, “Shearwave dispersion ultrasound vibrometry (SDUV) on swine kidney,” IEEE Transactions on Ultrasonics,
Ferroelectrics and Frequency Control, vol. 58, no. 12, 2011.
[122] F. G. Mitri, M. W. Urban, M. Fatemi, and J. F. Greenleaf, “Shear wave dispersion
ultrasonic vibrometry for measuring prostate shear stiffness and viscosity: an in vitro
pilot study,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 2, pp. 235–
242, 2011.
[123] C. Pislaru, M. W. Urban, I. Nenadic, and J. F. Greenleaf, “Shearwave dispersion ultrasound vibrometry applied to in vivo myocardium,” in Annual International Conference
of the IEEE Engineering in Medicine and Biology Society, pp. 2891–2894, IEEE, 2009.
[124] X. Zhang, R. R. Kinnick, M. Fatemi, and J. F. Greenleaf, “Noninvasive method for
estimation of complex elastic modulus of arterial vessels,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 52, no. 4, pp. 642–652, 2005.
[125] R. J. McGough, “Rapid calculations of time-harmonic nearfield pressures produced by
rectangular pistons,” The Journal of the Acoustical Society of America, vol. 115, no. 5,
pp. 1934–1941, 2004.
[126] J. A. Jensen, “Field: A program for simulating ultrasound systems,” in 10th
Nordicbaltic Conference on Biomedical Imaging, vol. 4, pp. 351–353, 1996.
[127] S. Holm, “Ultrasim-a toolbox for ultrasound field simulation,” University of Oslo, 2001.

131

[128] “http://www.signal.uu.se/toolbox/dream/userman.shtml,”
[129] J. A. Jensen and N. B. Svendsen, “Calculation of pressure fields from arbitrarily
shaped, apodized, and excited ultrasound transducers,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 39, no. 2, pp. 262–267, 1992.
[130] R. J. McGough, “Rapid calculations of time-harmonic nearfield pressures produced
by rectangular pistons,” Journal of the Acoustical Society of America, vol. 115, no. 5,
pp. 1934–1941, 2004.
[131] R. J. McGough, T. V. Samulski, and J. F. Kelly, “An efficient grid sectoring method
for calculations of the near-field pressure generated by a circular piston,” The Journal
of the Acoustical Society of America, vol. 115, no. 5, pp. 1942–1954, 2004.
[132] X. Zeng and R. J. McGough, “Evaluation of the angular spectrum approach for simulations of near-field pressures,” The Journal of the Acoustical Society of America,
vol. 123, no. 1, pp. 68–76, 2008.
[133] X. Zeng and R. J. McGough, “Optimal simulations of ultrasonic fields produced by
large thermal therapy arrays using the angular spectrum approach,” The Journal of
the Acoustical Society of America, vol. 125, no. 5, pp. 2967–2977, 2009.
[134] J. F. Kelly and R. J. McGough, “A time-space decomposition method for calculating
the nearfield pressure generated by a pulsed circular piston,” IEEE Transactions on
Ultrasonics, Ferroelectrics, and Frequency Control, vol. 53, no. 6, pp. 1150–1159, 2006.
[135] R. McGough and M. W. Urban, “Calculations of intensities for radiation force modeling
with the software package FOCUS,” in Proceedings of Meetings on Acoustics, vol. 9,
p. 020006, Acoustical Society of America, 2012.
[136] E. J. Alles, Y. Zhu, K. W. A. Van Dongen, and R. J. McGough, “Rapid transient
pressure field computations in the nearfield of circular transducers using frequencydomain time-space decomposition,” Ultrasonic Imaging, vol. 34, no. 4, pp. 237–260,
2012.
[137] K. R. Nightingale, R. W. Nightingale, M. L. Palmeri, and G. E. Trahey, “A finite
element model of remote palpation of breast lesions using radiation force: Factors
affecting tissue displacement,” Ultrasonic Imaging, vol. 22, no. 1, pp. 35–54, 2000.
[138] M. L. Palmeri, A. C. Sharma, R. R. Bouchard, R. W. Nightingale, and K. R. Nightingale, “A finite-element method model of soft tissue response to impulsive acoustic radiation force,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control,
vol. 52, no. 10, pp. 1699–1712, 2005.
[139] B. J. Fahey, K. R. Nightingale, R. C. Nelson, M. L. Palmeri, and G. E. Trahey,
“Acoustic radiation force impulse imaging of the abdomen: demonstration of feasibility
and utility,” Ultrasound in Medicine & Biology, vol. 31, no. 9, pp. 1185–1198, 2005.

132

[140] K. H. Lee, B. A. Szajewski, Z. Hah, K. J. Parker, and A. M. Maniatty, “Modeling
shear waves through a viscoelastic medium induced by acoustic radiation force,” International Journal for Numerical Methods in Biomedical Engineering, vol. 28, no. 6-7,
pp. 678–696, 2012.
[141] M. L. Palmeri, S. A. McAleavey, G. E. Trahey, and K. R. Nightingale, “Ultrasonic
tracking of acoustic radiation force-induced displacements in homogeneous media,”
IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 53, no. 7,
pp. 1300–1313, 2006.
[142] N. C. Rouze, M. H. Wang, M. L. Palmeri, and K. R. Nightingale, “Finite element
modeling of impulsive excitation and shear wave propagation in an incompressible,
transversely isotropic medium,” Journal of Biomechanics, vol. 46, no. 16, pp. 2761–
2768, 2013.
[143] M. Orescanin, Y. Wang, and M. F. Insana, “3-D FDTD simulation of shear waves for
evaluation of complex modulus imaging,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 58, no. 2, pp. 389–398, 2011.
[144] J. Langdon, K. Mercado, D. Dalecki, and S. McAleavey, “Compensating for Scholte
waves in single track location shearwave elasticity imaging,” The Journal of the Acoustical Society of America, vol. 137, no. 4, pp. 2364–2364, 2015.
[145] Y. Yang, M. Urban, and R. J. McGough, “A two-dimensional finite difference model
of shear wave propagation in anisotropic soft tissue,” in 2014 IEEE International
Ultrasonics Symposium, pp. 2323–2326, IEEE, 2014.
[146] B. Qiang, J. C. Brigham, S. Aristizabal, J. F. Greenleaf, X. Zhang, and M. W. Urban,
“Modeling transversely isotropic, viscoelastic, incompressible tissue-like materials with
application in ultrasound shear wave elastography,” Physics in Medicine and Biology,
vol. 60, no. 3, p. 1289, 2015.
[147] S. Callé, J.-P. Remenieras, O. B. Matar, M. E. Hachemi, and F. Patat, “Temporal
analysis of tissue displacement induced by a transient ultrasound radiation force,” The
Journal of the Acoustical Society of America, vol. 118, no. 5, pp. 2829–2840, 2005.
[148] 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.
[149] T. Deffieux, G. Montaldo, M. Tanter, and M. Fink, “Shear wave spectroscopy for
in vivo quantification of human soft tissues visco-elasticity,” IEEE Transactions on
Medical Imaging, vol. 28, no. 3, pp. 313–322, 2009.
[150] M. Couade, M. Pernot, C. Prada, E. Messas, J. Emmerich, P. Bruneval, A. Criton,
M. Fink, and M. Tanter, “Quantitative assessment of arterial wall biomechanical properties using shear wave imaging,” Ultrasound in Medicine & Biology, vol. 36, no. 10,
pp. 1662–1676, 2010.
133

[151] M. W. Urban, S. Chen, and M. Fatemi, “A review of shearwave dispersion ultrasound
vibrometry (SDUV) and its applications,” Current Medical Imaging Reviews, vol. 8,
no. 1, pp. 27–36, 2012.
[152] K. Aki and P. G. Richards, Quantitative Seismology. University Science Books, 2002.
[153] M. H. Sadd, Elasticity: Theory, Applications, and Numerics. Academic Press, 2009.
[154] S. Catheline, J.-L. Gennisson, G. Delon, M. Fink, R. Sinkus, S. Abouelkaram, and
J. Culioli, “Measurement of viscoelastic properties of homogeneous soft solid using
transient elastography: an inverse problem approach,” The Journal of the Acoustical
Society of America, vol. 116, no. 6, pp. 3734–3741, 2004.
[155] J. L. Rose, “A baseline and vision of ultrasonic guided wave inspection potential,”
Journal of Pressure Vessel Technology, vol. 124, no. 3, pp. 273–282, 2002.
[156] T. Loupas, J. T. Powers, and R. W. Gill, “An axial velocity estimator for ultrasound
blood flow imaging, based on a full evaluation of the doppler equation by means of a
two-dimensional autocorrelation approach,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 42, no. 4, pp. 672–688, 1995.
[157] H. Zhao, P. Song, M. W. Urban, J. F. Greenleaf, and S. Chen, “Shear wave speed
measurement using an unfocused ultrasound beam,” Ultrasound in Medicine & Biology,
vol. 38, no. 9, pp. 1646–1655, 2012.
[158] L. Sandrin, B. Fourquet, J. Hasquenoph, S. Yon, C. Fournier, F. Mal, C. Christidis,
M. Ziol, B. Poulet, F. Kazemi, M. Beaugrand, and R. Palau, “Transient elastography:
a new noninvasive method for assessment of hepatic fibrosis,” Ultrasound in Medicine
& Biology, vol. 29, no. 12, pp. 1705–1713, 2003.
[159] J. McLaughlin and D. Renzi, “Shear wave speed recovery in transient elastography
and supersonic imaging using propagating fronts,” Inverse Problems, vol. 22, no. 2,
p. 681, 2006.
[160] M. H. Wang, M. L. Palmeri, V. M. Rotemberg, N. C. Rouze, and K. R. Nightingale, “Improving the robustness of time-of-flight based shear wave speed reconstruction
methods using RANSAC in human liver in vivo,” Ultrasound in Medicine & Biology,
vol. 36, no. 5, pp. 802–813, 2010.
[161] N. C. Rouze, M. H. Wang, M. L. Palmeri, and K. R. Nightingale, “Robust estimation of
time-of-flight shear wave speed using a radon sum transformation,” IEEE Transactions
on Ultrasonics, Ferroelectrics and Frequency Control, vol. 57, no. 12, pp. 2662–2670,
2010.
[162] M. E. Hachemi, S. Callé, and J. Remenieras, “Transient displacement induced in shear
wave elastography: comparison between analytical results and ultrasound measurements,” Ultrasonics, vol. 44, pp. e221–e225, 2006.

134

[163] J. Ormachea, R. J. Lavarello, S. A. McAleavey, K. J. Parker, and B. Castaneda, “Shear
wave speed measurements using crawling wave sonoelastography and single tracking
location shear wave elasticity imaging for tissue characterization,” IEEE Transactions
on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 63, no. 9, pp. 1351–1360,
2016.
[164] R. S. Khalitov, S. N. Gurbatov, and I. Y. Demin, “The use of the verasonics ultrasound system to measure shear wave velocities in CIRS phantoms,” Physics of Wave
Phenomena, vol. 24, no. 1, pp. 73–76, 2016.
[165] Y. Yang, M. Urban, and R. J. McGough, “GPU-based green’s function simulations of
shear waves generated by an applied acoustic radiation force in elastic and viscoelastic
media,” Physics in Medicine & Biology, in review.
[166] H. Zhao, P. Song, M. W. Urban, R. R. Kinnick, M. Yin, J. F. Greenleaf, and S. Chen,
“Bias observed in time-of-flight shear wave speed measurements using radiation force
of a focused ultrasound beam,” Ultrasound in Medicine & Biology, vol. 37, no. 11,
pp. 1884–1892, 2011.
[167] I. Z. Nenadic, M. W. Urban, H. Zhao, W. Sanchez, P. E. Morgan, J. F. Greenleaf, and
S. Chen, “Application of attenuation measuring ultrasound shearwave elastography in
8 post-transplant liver patients,” in 2014 IEEE International Ultrasonics Symposium,
pp. 987–990, IEEE, 2014.
[168] M. Couade, M. Pernot, E. Messas, A. Bel, M. Ba, A. Hagege, M. Fink, and M. Tanter, “In vivo quantitative mapping of myocardial stiffening and transmural anisotropy
during the cardiac cycle,” IEEE Transactions on Medical Imaging, vol. 30, no. 2,
pp. 295–305, 2011.
[169] C. Hazard, Z. Hah, D. Rubens, and K. Parker, “Integration of crawling waves in an
ultrasound imaging system. part 1: system and design considerations,” Ultrasound in
Medicine & Biology, vol. 38, no. 2, pp. 296–311, 2012.
[170] Z. Hah, C. Hazard, B. Mills, C. Barry, D. Rubens, and K. Parker, “Integration of
crawling waves in an ultrasound imaging system. part 2: signal processing and applications,” Ultrasound in Medicine & Biology, vol. 38, no. 2, pp. 312–323, 2012.
[171] M. L. Palmeri and K. R. Nightingale, “On the thermal effects associated with radiation
force imaging of soft tissue,” IEEE Transactions on Ultrasonics, Ferroelectrics and
Frequency Control, vol. 51, no. 5, pp. 551–565, 2004.
[172] F. Varray, C. Cachard, A. Ramalli, P. Tortoli, and O. Basset, “Simulation of ultrasound nonlinear propagation on GPU using a generalized angular spectrum method,”
EURASIP Journal on Image and Video Processing, vol. 2011, no. 1, pp. 1–6, 2011.
[173] B. E. Treeby, M. Tumen, and B. T. Cox, “Time domain simulation of harmonic ultrasound images and beam patterns in 3D using the k-space pseudospectral method,”
in 2011 Medical Image Computing and Computer-Assisted Intervention, pp. 363–370,
Springer, 2011.
135

[174] P. M. Novotny, J. A. Stoll, N. V. Vasilyev, J. Pedro, P. E. Dupont, T. E. Zickler,
and R. D. Howe, “GPU based real-time instrument tracking with three-dimensional
ultrasound,” Medical Image Analysis, vol. 11, no. 5, pp. 458–464, 2007.
[175] S. Rosenzweig, M. Palmeri, and K. Nightingale, “GPU-based real-time small displacement estimation with ultrasound,” IEEE Transactions on Ultrasonics, Ferroelectrics,
and Frequency Control, vol. 58, no. 2, pp. 399–405, 2011.
[176] J. P. Asen, J. I. Buskenes, C. I. C. Nilsen, A. Austeng, and S. Holm, “Implementing
capon beamforming on a GPU for real-time cardiac ultrasound imaging,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 61, no. 1, pp. 76–85,
2014.
[177] J. Chen, B. Yiu, H. K. H. So, and A. C. H. Yu, “Real-time GPU-based adaptive beamformer for high quality ultrasound imaging,” in 2011 IEEE International Ultrasonics
Symposium, pp. 474–477, IEEE, 2011.
[178] P. Song, A. Manduca, H. Zhao, M. W. Urban, J. F. Greenleaf, and S. Chen, “Fast
shear compounding using robust 2-d shear wave speed calculation and multi-directional
filtering,” Ultrasound in Medicine & Biology, vol. 40, no. 6, pp. 1343–1355, 2014.
[179] S. L. Friedman, “Liver fibrosis–from bench to bedside,” Journal of Hepatology, vol. 38,
pp. 38–53, 2003.
[180] M. J. Buckingham, “Causality, stokes wave equation, and acoustic pulse propagation
in a viscous fluid,” Physical Review E, vol. 72, no. 2, p. 026610, 2005.
[181] D. Klatt, U. Hamhaber, P. Asbach, J. Braun, and I. Sack, “Noninvasive assessment
of the rheological behavior of human organs using multifrequency MR elastography:
a study of brain and liver viscoelasticity,” Physics in Medicine and Biology, vol. 52,
no. 24, p. 7281, 2007.
[182] I. Sack, K. Jöhrens, J. Würfel, and J. Braun, “Structure-sensitive elastography: on
the viscoelastic powerlaw behavior of in vivo human tissue in health and disease,” Soft
Matter, vol. 9, no. 24, pp. 5672–5680, 2013.
[183] C. Coussot, S. Kalyanam, R. Yapp, and M. F. Insana, “Fractional derivative models for ultrasonic characterization of polymer and breast tissue viscoelasticity,” IEEE
Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 56, no. 4,
pp. 715–726, 2009.
[184] L. S. Taylor, A. L. Lerner, D. J. Rubens, and K. J. Parker, “A kelvin-voight fractional
derivative model for viscoelastic characterization of liver tissue,” in ASME 2002 International Mechanical Engineering Congress and Exposition, pp. 447–448, American
Society of Mechanical Engineers, 2002.
[185] 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 & Biology, vol. 33, no. 10,
pp. 1617–1631, 2007.
136

[186] K. Hoyt, B. Castaneda, M. Zhang, P. Nigwekar, P. A. di Sant’Agnese, J. V. Joseph,
J. Strang, D. J. Rubens, and K. J. Parker, “Tissue elasticity properties as biomarkers
for prostate cancer,” Cancer Biomarkers, vol. 4, no. 4, 5, pp. 213–225, 2008.
[187] M. Zhang, P. Nigwekar, B. Castaneda, K. Hoyt, J. V. Joseph, A. di Sant’Agnese,
J. G. Messing, E. M .and Strang, D. J. Rubens, and K. J. Parker, “Quantitative
characterization of viscoelastic properties of human prostate correlated with histology,”
Ultrasound in Medicine & Biology, vol. 34, no. 7, pp. 1033–1042, 2008.
[188] D. Craiem, F. J. Rojo P., J. M. Atienza R., G. V. Guinea, and R. L. Armentano,
“Fractional calculus applied to model arterial viscoelasticity,” Latin American Applied
Research, vol. 38, no. 2, pp. 141–145, 2008.

137