Comparative study of shear wave-based elastography techniques in optical coherence tomography

Fernando Zvietcovich; Jannick P. Rolland; Jianing Yao; Panomsak Meemon; Kevin J. Parker

doi:10.1117/1.JBO.22.3.035010

30 March 2017 Comparative study of shear wave-based elastography techniques in optical coherence tomography

Fernando Zvietcovich, Jannick P. Rolland, Jianing Yao, Panomsak Meemon, Kevin J. Parker

Author Affiliations +

Journal of Biomedical Optics, Vol. 22, Issue 3, 035010 (March 2017). https://doi.org/10.1117/1.JBO.22.3.035010

Abstract

We compare five optical coherence elastography techniques able to estimate the shear speed of waves generated by one and two sources of excitation. The first two techniques make use of one piezoelectric actuator in order to produce a continuous shear wave propagation or a tone-burst propagation (TBP) of 400 Hz over a gelatin tissue-mimicking phantom. The remaining techniques utilize a second actuator located on the opposite side of the region of interest in order to create three types of interference patterns: crawling waves, swept crawling waves, and standing waves, depending on the selection of the frequency difference between the two actuators. We evaluated accuracy, contrast to noise ratio, resolution, and acquisition time for each technique during experiments. Numerical simulations were also performed in order to support the experimental findings. Results suggest that in the presence of strong internal reflections, single source methods are more accurate and less variable when compared to the two-actuator methods. In particular, TBP reports the best performance with an accuracy error <4.1%. Finally, the TBP was tested in a fresh chicken tibialis anterior muscle with a localized thermally ablated lesion in order to evaluate its performance in biological tissue.

1. Introduction

Elastic properties of tissue provide valuable information regarding pathological processes such as liver fibrosis,¹ arteriosclerosis,²^,³ breast cancer,⁴ and glaucoma,⁵ according to a number of studies reported during the last three decades. For instance, ultrasound elastography (USE) is used to measure the biomechanical properties of tissue with a resolution and penetration depth on the order of 1, and 10 mm, respectively.⁶ Magnetic resonance elastography, on the other hand, has a wider scanning range ( $\sim 10$ to 100 mm) with resolutions on the order of $\sim 10 mm$ .⁷ Optical coherence tomography (OCT) is an imaging modality used to acquire structural images of tissue ranging from 1 to 10 mm in the lateral dimension with micrometric resolution.⁸ The posterior application of elastography in OCT, commonly known as optical coherence elastography (OCE), introduced a new level of tissue analysis due to its unique capability in resolution (microscale spatial resolution), which enables the elastography imaging of small scale biocomponents.⁹

In general, the measurement of the elastic properties of tissue using OCE requires the application of force to the sample in order to generate localized tissue displacements. The level of tissue displacement will give information about the elastic (and viscous) properties of the sample. Depending on the mechanical loading characteristics, OCE techniques can be classified into two categories: quasi-static (compression), and dynamic (transient excitation or sinusoidal steady state).⁹ Quasi-static OCE was first demonstrated by Schmitt.¹⁰ It is also important to mention the work of Wang et al.¹¹ in which Doppler OCT is used for detecting the tissue motion produced by a translational stage in the sample.

Dynamic OCE enables the quantitative measurement of the biomechanical properties of tissue by estimating the speed of shear waves propagating in the sample without any previous knowledge of the loading stress.¹² Sinusoidal harmonic and impulse excitation have been applied in dynamic OCE, in which the perturbation wavefront propagates along the sample. Subsequently, the different local mechanical properties of the sample will modify the wave speed that can be measured and interpreted using an appropriate biomechanical model. Several investigations were conducted in this area when mechanical actuators,¹²^–¹⁵ acoustic radiation force,¹⁶^,¹⁷ air-puff excitation,¹⁸^–²⁰ and magnetic transducers²¹ were utilized in the wave generation. In particular, Manapuram et al.¹³ used a swept-source OCT (SS-OCT) for measuring shear wave pulse speed in gelatin phantoms. Similarly, Song et al.¹⁴ proved that a phase-sensitive OCT system can be used for tracking mechanical waves in gelatin samples at a frame rate of 47 kHz created with a piezoelectric transducer. Li et al.¹⁵ used the vibration method and surface acoustic waves (SAW) to determine the biomechanical properties of epidermal layers of skin. Nahas et al.¹⁶ induced shear waves along the sample using acoustic radiation force and measured its phase speed in order to generate elasticity maps in tissue-mimicking phantoms and ex vivo rat brain. Air-puff excitation was utilized for generating lamb waves (a type of SAW) in ex vivo rabbit cornea by Wang and Larin¹⁸ in order to distinguish four biomechanically distinct layers. Finally, the measurement of viscoelasticity in elastic and viscous phantoms as well as in biological tissues by calculating the shear wave speed produced by localized inclusion of magnetic nanoparticles was proven by Ahmad et al.²¹

Given the many approaches, we inquire about what combination of sources and analyses will produce the most accurate, rapid, and high-resolution images of tissue stiffness. In this study, we focus on shear wave techniques associated with dynamic OCE, with both transient tone-burst and sinusoidal steady-state mechanical excitation generated by one or two external sources. Some leading estimators are applied to the displacement data to produce local estimates of shear wave speed, which are related to the stiffness of the material and are the basis of elastography images. In the remainder of this paper, we will summarize in Sec. 2 the theory of shear wave propagation (SWP) and measurements using phase-sensitive OCT. In Sec. 3, we will detail experiments conducted with the various techniques, and results that compare findings will be presented in Sec. 4. Finally, we will discuss the results and conclude this paper in Sec. 5.

2. Theory

This section describes the principles of how motion is measured using phase-sensitive OCT, and the physics related to the propagation of shear waves in a linear elastic medium generated by one or two external sources, including a special case of SAW. Phase velocity estimation of waves using a phase-sensitive OCT system will be also discussed.

2.1.

Phase-Sensitive Optical Coherence Tomography

Phase-sensitive OCT (PhS-OCT) measures local minute displacement of particles based on the phase difference between two consecutive complex OCT signals reconstructed from Fourier transforms (FT) of the detected A-line OCT interference spectra that are collected continuously at approximately the same lateral position on the sample. In this research, the phase sensitive detection was implemented on an SS-OCT system with a Mach–Zehnder interferometer configuration [see Fig. 1(a)]. Detailed specifications of the system will be provided in Sec. 3.2. Methodologically, in a medium with refractive index $n$ , the phase difference $Δ ϕ (z) = ϕ (z, t_{1}) - ϕ (z, t_{0})$ at two consecutive instants $t_{0}$ and $t_{1}$ ( $t_{0} < t_{1}$ ) for a given lateral position $x_{0}$ is related to the particle displacement in the axial direction by

Eq. (1)

u_{z} (z) = Δ ϕ (z) \frac{λ_{0}}{4 π n},

where

ϕ

is the phase of the acquired A-line signal,

n

is the refractive index of the medium, and

λ_{0}

is the center wavelength of the laser. Figure 1(b) shows the OCT system, in which the OCT beam penetrates the sample in the axial direction in order to calculate the displacement along the A-line signal. The repetition of this procedure for a range of lateral positions within a region of interest (ROI) yields a two-dimensional (2-D) displacement map

u_{z} (x, z)

, in which a mechanical wave (i.e., shear wave) propagating in the medium can be tracked.

Fig. 1

(a) Sweep-source OCT system. (b) OCT laser beam penetration in sample.

PhS-OCT possesses a high displacement sensitivity, on the order of nanometers,¹¹ with respect to other motion detection techniques such as speckle tracking. However, PhS-OCT may be subjected to speckle decorrelation, motion artifacts, and noise. These negative effects can be attenuated by using a spatial averaging filter. In addition, artifacts in $u_{z} (z)$ created by the motion of the air–medium interface can be compensated as discussed in the work of Song et al.¹⁴

2.1.1.

Motion detection

The estimation of $Δ ϕ (z)$ yields the direct calculation of displacement $u (z)$ in the sample. There exist various methods for the estimation of $Δ ϕ (z)$ , also called Doppler methods, which vary in accuracy and noise reduction. Consider two consecutive complex-valued A-line signals

Eq. (2)

I (z, t_{0}) = | I (z, t_{0}) | e^{j ϕ (z, t_{0})},

Eq. (3)

I (z, t_{1}) = | I (z, t_{1}) | e^{j ϕ (z, t_{1})},

where

| I |

represents the magnitude,

ϕ

is the phase of the signal at two consecutive instants

t_{0}

and

t_{1}

, and

t_{0} < t_{1}

. A direct way to estimate the phase difference

Δ ϕ (z) = ϕ (z, t_{1}) - ϕ (z, t_{0})

can be achieved by

Eq. (4)

Δ ϕ (z) = \arg [I (z, t_{1}) I {(z, t_{0})}^{*}],

Eq. (5)

Δ ϕ (z) = \arctan {\frac{Im [I (z, t_{1}) I {(z, t_{0})}^{*}]}{Re [I (z, t_{1}) I {(z, t_{0})}^{*}]}},

where

I^{*}

is the complex conjugated version of

I

, and

Im {\cdot}

and

Re {\cdot}

are the imaginary and real-part operators of a given vector, respectively. The main problem with this approach lies in its high instability and sensitivity to noise. Based on the study developed by Pinton et al.,²² we chose the Loupas et al.²³ 2-D autocorrelation as the motion estimation method. Loupas’ algorithm uses depth and lateral samples to calculate displacement. When we take two samples within the axial range,

Δ ϕ (z)

can be represented as

Eq. (6)

Δ ϕ (z) = \frac{\arctan {\frac{\sum_{m = 0}^{M - 1} [Q (z_{m}, t_{0}) I I (z_{m}, t_{1}) - I I (z_{m}, t_{0}) Q (z_{m}, t_{1})]}{\sum_{m = 0}^{M - 1} [I I (z_{m}, t_{0}) I I (z_{m}, t_{1}) + Q (z_{m}, t_{0}) Q (z_{m}, t_{1})]}}}{1 + \arctan {\frac{\sum_{m = 0}^{M - 2} \sum_{n = 0}^{1} [Q (z_{m}, t_{n}) I I (z_{m + 1}, t_{n}) - I I (z_{m}, t_{n}) Q (z_{m + 1}, t_{n})]}{\sum_{m = 0}^{M - 2} \sum_{n = 0}^{1} [I I (z_{m}, t_{n}) I I (z_{m + 1}, t_{n}) + Q (z_{m}, t_{n}) Q (z_{m + 1}, t_{n})]}} / 2 π},

where

z_{m}

is the element

m

’th along the

z

-axis,

I I (z, t) = Re {I (z, t)}

, and

Q (z, t) = Im {I (z, t)}

, and

M

is the size of the window in the axial direction in which the phase estimation is performed for a given

z

. For an intensity vector

I

of 1025 elements, a window size

M = 20

significantly reduces the noise without a strong resolution loss. Applying Eq. (6), the displacement can be calculated more accurately than with the classic method [Eq. (5)].

2.2.

Shear Waves in Elastic Media

PhS-OCT and motion detection methods enable the calculation of 2-D displacement frames describing wave propagation in a sample as described in Sec. 2.1. Now, the physical relationship between shear wave speed and shear modulus needs to be established for two cases: one-source and two-sources of vibration.

2.2.1.

One-source vibration elastography

A shear wave is a transverse body wave whose displacement is perpendicular to the direction of the wave propagation. The propagation of a shear and harmonic plane wave modeled in the $x z$ -plane [as described in Fig. 1(b)] and traveling toward the $x$ -axis in a linear elastic medium is governed by the one-dimensional (1-D) solution of the wave equation given by

Eq. (7)

W (x, z; t) = A e^{- α (x + \frac{D}{2})} e^{- j [k (x + \frac{D}{2}) - ω t]},

where

W

is defined as the displacement along the

z

-axis,

A

is the amplitude of the wave,

α

is the attenuation coefficient of the medium,

k

is the wave number,

ω

is the angular frequency, and

x = - D / 2

is the location where the wave is generated.²⁴ Then, at a given time

t_{0}

, the particle displacement generated by the shear wave [Fig. 2(a)] can be measured by a PhS-OCT system using Eq. (1), and

u_{z} (x, z) = W (x, z; t_{0})

. The phase velocity of the shear wave, also referred to as the speed of the shear wave, is given by

Eq. (8)

c_{shear} = \frac{ω}{k} = λ f,

where

λ

is the wavelength of the mechanical excitation and

f

is the frequency of the shear wave. The shear wave speed is related to the shear modulus

G

and mass density

ρ

of the medium as⁶

Eq. (9)

c_{shear} = \sqrt{\frac{G}{ρ}} .

Fig. 2

(a) SWP in a solid elastic medium [Eq. (7)]. (b) Interference of two waves traveling in opposite directions [Eqs. (12) and (13)]. (c) CrW interference pattern and its amplitude profile [Eq. (14)].

In general, biological tissues are nearly incompressible and the Poisson’s ratio can be approximated to $ν = 0.5$ .⁶ Then the shear modulus $G$ can be related to the elastic modulus (Young’s modulus) $E$ as $E = 3 G$ .⁶ Therefore, the Young’s modulus $E$ can be determined by measuring the speed of the shear wave in a homogenous isotropic medium. In dynamic OCE applications, the frequency of the excitation is known; then the calculation of $c_{shear}$ is reduced to the estimation of $λ$ or, equivalently, by estimating the wave number $k = 2 π / λ$ . In addition, the shear wave speed can also be estimated by using the kinematic relationship

Eq. (10)

c_{shear} = \frac{Δ d_{shear}}{Δ t_{shear}},

where

Δ d_{shear}

is the propagated distance of a transient shear wave tone during time

Δ t_{shear}

. This estimation technique is also called the time-of-flight method and is widely used in ultrasound (US) and OCT elastography.²⁵ Shear wave velocity in tissue typically ranges from 0.1 to

10 m / s

;²⁶ therefore, OCT systems with a lateral field of view on the order of 30 mm are not able to fully track a single shear wavefront due to the need of scanning line-by-line with repetitions for Doppler motion estimates. Several approaches solving this estimation problem are proposed and presented in Sec. 3.

In addition, SAW are produced when a vibration is generated at a solid-vacuum or solid–fluid interface.²⁷ Since the penetration depth of an OCT system typically ranges some millimeters from the interface, SAW are more likely to be scanned than shear waves. If we consider a semi-infinite solid medium, the predominant SAW propagating in the solid–vacuum interface is Rayleigh waves. The relationship between shear wave and Rayleigh wave phase speed in a linear isotropic medium for a Poisson’s ratio $ν = 0.5$ is approximately given by²⁷

Eq. (11)

c_{shear} \approx 1.05 * c_{Rayleigh} .

For further details about Rayleigh waves in OCE, we refer the reader to our previous study by Zvietcovich et al.²⁷

2.2.2.

Two-source vibration elastography

When two vibration sources produce shear waves traveling in opposite directions [Fig. 2(b)], an interference pattern with interesting properties is created. Let us assume two plane waves propagating in the right and left directions $W_{left}$ and $W_{right}$ parallel to the $x$ -axis, and with different frequencies $ω$ and $ω + Δ ω$ , respectively. For a selected depth $z_{0}$ , we can describe the propagation of those waves in the $x - t$ domain as

Eq. (12)

W_{left} (x, t) = A e^{- α (x + \frac{D}{2})} e^{- j [k (x + \frac{D}{2}) - ω t + φ_{l}]},

Eq. (13)

W_{right} (x, t) = A e^{- α (\frac{D}{2} - x)} e^{- j [(k + Δ k) (\frac{D}{2} - x) - (ω + Δ ω) t + φ_{r}]},

where

k

and

k + Δ k

are the wave numbers of

W_{left}

and

W_{right}

, respectively;

φ_{l}

and

φ_{r}

are the initial arbitrary phases, and

D

is the separation distance between the two sources that are assumed to be located at

x = \pm D / 2

, respectively. According to Wu et al.,²⁸ assuming no reflections in the medium, the square amplitude of the interference

U (x, t) = W_{left} + W_{right}

is given by

Eq. (14)

{| U (x, t) |}^{2} = 2 A^{2} e^{- α D} {\cosh (2 α x) + \cos [Δ ω t + 2 (k + \frac{Δ k}{2}) x - φ_{0}]} .

The hyperbolic cosine in Eq. (14) refers to the attenuation of the medium, and the cosine term refers to a slow speed wave, also called a crawling wave (CrW) pattern²⁸ [Fig. 2(c)]. The phase term

φ_{0} = \frac{Δ k}{2} D + Δ φ

is constant and

Δ φ = φ_{r} - φ_{l}

. The CrW speed

c_{Cr W}

is related to the shear wave speed

c_{shear}

by

Eq. (15)

c_{CrW} = c_{shear} \frac{Δ ω}{2 ω + Δ ω},

where

c_{shear} = ω / k

, and

c_{Cr W} = Δ ω / (2 k + Δ k)

. In US and OCT vibration elastography applications, the parameter

Δ ω

is selected according to the hardware acquisition velocity requirements. There exist various methods for estimating

c_{shear}

out of the CrW pattern and these will be discussed in Sec. 3. Of particular interest is the case in which both vibration sources have the same angular frequency

ω

, that is,

Δ ω = 0

. Then,

c_{CrW} = 0

, and the CrW pattern becomes a standing wave (SW). In this situation,

Eq. (16)

{| U (x, t) |}^{2} = 2 A^{2} e^{- α D} [\cosh (2 α x) + \cos (2 k x - φ_{0})],

and the shear wave speed estimation reduces to the estimation of

k

in the SW pattern.

2.3.

Shear Wave Phase-Velocity Estimators

As discussed in Sec. 2.2, the localized phase velocity of a plane shear wave propagating in an elastic medium reduces to the calculation of the time-of-flight of a transient tone-burst produced by one vibration source, or the spatially dependent wave number $k (x)$ of a harmonic steady-state vibration produced by one or two sources. Based on these approaches, three common techniques to estimate $c_{shear}$ from the displacement field measured by PhS-OCT will be presented and discussed in this section.

2.3.1.

Phase derivative estimator

This estimator allows the shear wave speed to be determined based on the calculation of space–time representations. This is done by calculating the instantaneous phase of each harmonic signal at a given lateral and depth position.²⁹ For a given depth $z_{0}$ , Eq. (14) is a harmonic space–time dependent 2-D signal. Using a finite impulse response filter on the time-domain, we extract the component of the signal centered at the angular frequency $Δ ω$ in order to remove the nonharmonic cosh term. Then, assuming $Δ ω ≪ ω$ , it is also true that $Δ k ≪ k$ , and

Eq. (17)

{| U (x, t) |}^{2} = B {\cos [Δ ω t + ϕ (x)]},

where

B = 2 A^{2} e^{- α D}

is the amplitude and

ϕ (x) = 2 k x - φ_{0}

is the phase signal for each lateral position. The estimated version of the phase

\tilde{ϕ} (x)

can be computed by means of harmonic fitting or by using the FT approach. The estimated shear wave speed

{\tilde{c}}_{shear} (x)

can be calculated using

Eq. (18)

{\tilde{c}}_{shear} (x) = \frac{ω}{\tilde{k}} = \frac{4 π f}{{\tilde{ϕ}}^{'} (x)},

where

{\tilde{ϕ}}^{'} (x)

is the spatial derivative of

\tilde{ϕ} (x)

, and

ω = 2 π f

. The derivative operator is vulnerable to noise corruption; then a smoothing procedure can be applied to

\tilde{ϕ} (x)

at the expense of losing lateral resolution.

{\tilde{c}}_{shear} (x)

is computed for each depth

z

in order to create a 2-D shear wave velocity map.

2.3.2.

Hoyt’s estimator

The Hoyt’s estimator is able to calculate the spatial frequency of a wave based on a complex cross-correlation procedure.³⁰ The signal of Eq. (17) is converted to a complex-valued analytic signal $s_{A}$ of $N$ elements using Hilbert transform. Subsequently, using a kernel window of size $L$ , the first lag of the autocorrelation function $R$ is

Eq. (19)

R (n, 1) = \frac{1}{L - 1} \sum_{m = 0}^{L - 2} s_{A}^{*} (n + m) s_{A} (n + m + 1) .

Therefore, the spatial wave number can be estimated as

Eq. (20)

k (x_{n}) = \frac{1}{Δ_{x}} \arctan {\frac{Im [R (n, 1)]}{Re [R (n, 1)]}},

where

x_{n} = n Δ_{x}

, for

n = 0, 1, 2, \dots, N - 1 - L

and

Δ_{x}

is the lateral sampling resolution. Using Eqs. (20) and (18), a 2-D shear-wave-speed map (SWSM) can be computed.

2.3.3.

Peak tracking for time-of-flight estimator

This estimator is able to calculate the laterally resolved wave speed of a transient tone-burst propagating in an elastic medium by tracking the main peak of the burst in space and time. Herein, the local group velocity of the tone-burst is approximated by the slope of the space–time curve representing the peak trajectory within a window of size $w$ . Equation (10) can be adapted to a least squares regression approach in which lateral spatial and time samples within a window are stored in $p (n) = [x_{n}, x_{n + 1}, x_{n + 2}, \dots, x_{n + w - 1}]$ , and $q (n) = [t_{n}, t_{n + 1}, t_{n + 2}, \dots, t_{n + w - 1}]$ , respectively; and the laterally resolved group speed is estimated as

Eq. (21)

c_{peak} (x_{n}) = \frac{[\sum_{i = 0}^{w - 1} p_{i} (n) q_{i} (n)] - w \bar{p (n)} \bar{q (n)}}{[\sum_{i = 0}^{w - 1} q_{i} {(n)}^{2}] - w {\bar{q (n)}}^{2}},

where

\bar{q (n)} = \frac{1}{w} \sum_{i = 0}^{w - 1} q_{i} (n)

,

\bar{p (n)} = \frac{1}{w} \sum_{i = 0}^{w - 1} p_{i} (n)

,

x_{n} = n Δ_{x}

,

t_{n} = n Δ_{t}

for

n = 0, 1, 2, \dots, N - w

, and

Δ_{x}

,

Δ_{t}

are the lateral spatial and time sampling resolution, respectively. A larger size of

w

will increase the accuracy of the estimation at the cost of losing resolution in the lateral axis.

3. Materials and Methods

In this section, tissue-mimicking phantom preparation and measurement, the setup of the PhS-OCT system, and five dynamic OCE methods are described. These methods are based on the use of one or two vibration sources and will be discussed in terms of numerical simulations and experiments.

3.1.

Sample Preparation and Measurement

A two-sided gelatin phantom was prepared by using 10% and 15% gelatin powder concentration in each side in order to ensure differentiated mechanical properties (Young’s modulus). On both sides, the same concentration of 2% of intralipid powder (coffee creamer) was used in order to provide the phantom with light scattering properties. The refractive index of the phantom was $\sim 1.35$ , which is close to that of real tissue. The shape of the phantom is a rectangular cuboid of dimensions $15 \times 7 \times 3 {cm}^{3}$ [Fig. 3(a)].

The Young’s modulus of each side of the phantom was measured by taking three cylindrical samples ( $n = 3$ ) $\sim 45 - mm$ in diameter, and measured using a stress–relaxation compression test. The measurement was conducted by using an MTS Q-test/5 universal testing machine (MTS, Eden Prairie, Minnesota) with a 50 N load cell using a compression rate of $0.5 mm / s$ , a strain value of 5%, and total measurement time of 600 s. The stress–time plots obtained by the machine were fitted to a fractional derivative standard linear solid (FD-SLS) model in order to obtain a frequency-dependent Young’s modulus according to Schmidt and Gaul.³¹ Finally, using Eq. (9) for a material density of $1 g / {cm}^{3}$ , and assuming a nearly incompressible homogeneous and isotropic material (Poisson’s ratio of 0.5), the shear wave speed for each side of the phantom was estimated and used as ground truth measurement.

For the biological tissue test, a tibialis anterior muscle was dissected from a fresh roaster chicken and subjected to a localized thermal ablation in the surface in order to produce a stiffer lesion when compared to normal muscle. The muscle was not subjected to any external force in order to prevent passive muscle resistance effect.

3.2.

Phase-Sensitive Optical Coherence Tomography Experimental Setup

A PhS-OCT system was implemented utilizing a swept source (HSL-2100-WR, Santec, Aichi, Japan) with a center wavelength of 1318 nm and a full-width half-maximum (FWHM) bandwidth of 125 nm. The frequency sweep rate of the light source was 20 kHz. The source power that entered the OCT interferometer was split by a 10/90 fiber coupler into the reference and sample arms, respectively. In the reference arm, a custom Fourier domain optical delay line was used for dispersion compensation. In the sample arm, a collimated light beam diameter of 6.7 mm at $1 / e^{2}$ was directed onto a test phantom by a focusing imaging lens (LSM05, Thorlabs Inc., New Jersey), coupled with a galvanometer scanning mirror placed at the front focal plane of the imaging lens to achieve flat-field lateral scanning. The back-scattered light from the sample was recombined with the light reflected from the reference mirror with a 50/50 fiber coupler. The time-encoded spectral interference signal was detected by a balanced photodetector (1817-FC, New Focus, California), and then digitized with a $500 MSamples / s$ , 12-bit-resolution analog-to-digital converter (ATS9350, AlazarTech, Québec, Canada). The maximum sensitivity of the system was measured to be 112 dB.³² The imaging depth of the system was measured to be 5 mm in air ( $- 10 dB$ sensitive fall-off). The optical lateral resolution was $\sim 30 μ m$ , and the FWHM of the axial point spread function after dispersion compensation was $10 μ m$ . The synchronized control of the galvanometer and the OCT data acquisition was conducted through a LabView platform connected to a work station. The phase stability of the system was calculated as the standard deviation of the temporal fluctuations of the Doppler phase-shift ( $Δ ϕ_{err}$ ) while imaging a static structure.³³ Result shows $Δ ϕ_{err} = 4.6 mrad$ when using the Loupas’ algorithm of Eq. (6). The displacement sensitivity is measured as the minimum detectable axial displacement ( $u_{z, \min}$ ) by plugging $Δ ϕ_{err}$ in Eq. (1). We found a displacement of $u_{z, \min} = 0.358 nm$ . Finally, the maximum axial displacement supported by the system without unwrapping the phase-shift signal ( $Δ ϕ_{\max} = π$ ) is $u_{z, \max} = 0.24 μ m$ . Typical values of SNR of the OCT signal found in the gelatin phantoms were in the range of 19–23 dB. The SNR was obtained by calculating the standard deviation of the noise in a region of 0.2 × 1 mm right above the phantom surface which is located at the approximate depth of 0.25 mm in the ROI, and the average intensity of the phantom signal right below the phantom surface using the same region size.

The separation between the phantom and the surface of the objective lens was $\sim 9 cm$ . Two piezoelectric actuators (APC 40-2020, APC International, Ltd., Macheyville, Pennsylvania) were located on each side of the phantom as shown in Fig. 3(b). These actuators were connected to stereo amplifiers (LP-2020A+, Lepai, Houston, Texas) that receive harmonic signals from a dual channel function generator (AFG320, Tektronix). The maximum lateral field of view of the ROI in the phantom was 30 mm and the maximum depth was 2.5 mm. Both actuators were located $\sim 8 cm$ from the center of the ROI [Fig. 3(b)]. The out-of-plane vertical vibration of the actuators produced Rayleigh waves that were approximately planar within the ROI. The vibration of each actuator was adjusted to have the same amplitude in the ROI’s center in order to ensure maximum contrast, and to have a displacement magnitude $< 3 μ m$ in order to avoid nonlinear effects in elastic behavior of the material. In addition, the maximum displacement that can be detected by the OCT system is $u_{z, \max} = 0.24 μ m$ without unwrapping the Doppler phase signal. Since waves travel $\sim 6 cm$ from the excitation point to the ROI, waves will be attenuated by the medium so that wave displacement amplitudes within the ROI are in the full dynamic range of displacement sensitivity of the OCT system. The excitation frequency, acquisition protocol, and the use of one or two vibration sources will depend on the dynamic OCE method. For all cases, we use Eq. (11) to convert Rayleigh wave speed to shear wave speed.

Fig. 3

(a) Two-sided gelatin phantom (10%, 15% gel concentration each side). (b) Experimental setup showing location of vibrators, ROI, and dimensions.

3.3.

Dynamic Optical Coherence Elastography Methods

Five dynamic OCE methods are considered in this comparative study. Each method is classified by a vibration source configuration, OCT acquisition protocol, which will be described with each method, and the phase speed estimator as shown in Fig. 4. For all cases, the Loupas algorithm was used for the acquisition of motion frames, and the SAW speed calculated using phase speed estimators was corrected using the Rayleigh-shear speed relationship [Eq. (11)].

Fig. 4

Flowchart describing the classification of five dynamic OCE methods considered in the comparative study. (a) General description of the classification criteria of an OCE method. (b) Description of the algorithms and configurations of the five dynamic OCE methods corresponding to the classification criteria shown in the top chart.

Three differentiated OCT acquisition protocols were used in each method accordingly: M-mode, M–B mode, and B-raster mode. M-mode, also called motion mode, acquires a sequence of A-scans at every spatial location. If just one A-scan is acquired for each lateral position covering the entire ROI, we say that a B-mode image has been generated. M–B mode performs the M-mode scanning synchronized with the repeated firing of vibration waves in the sample. Then, at the end, the data can be reorganized as a collection of B-mode frames. Finally, the B-raster mode generates fast B-scan frames by continuously acquiring A-scans along the lateral axes until the ROI is covered. For each OCE method, the corresponding protocol will be described in detail.

Phase speed estimators are used according to the source configuration and acquisition protocol. Phase derivative estimator performs more accurately than Hoyt’s estimator for OCE methods using a sinusoidal steady-state vibration. Then we will preferably use a phase derivative estimator when possible; otherwise, the Hoyt’s method will be utilized. Peak tracking estimator is intended to estimate the time of flight of a transient tone-burst. Therefore, this estimator is reserved for methods using a single source producing a tone-burst.

3.3.1.

Crawling waves method

This method was initially proposed for US elastography in order to reduce the acquisition speed requirement of US scanners and still be able to measure the speed of shear waves in tissue.²⁸ The same idea was applied for the first time in OCE as reported by Meemon et al.³⁴ The two harmonic actuators are set at frequencies $f_{left} = 400 Hz$ and $f_{right} = 400.4 Hz$ , respectively, yielding a $Δ f = 0.4 Hz$ or $Δ ω \approx 2.51 rad / s$ . In the CrW method, the acquisition protocol is the B-raster mode that consists of scanning multiple B-mode motion frames at a frame rate of $5 frames / s$ using Loupas’ approach²³ described in Sec. 2.1.1. Each frame is formed by 1500 A-lines covering a lateral field of view of 10 mm with a lateral sampling resolution of $6.6 μ m$ , as described in Fig. 5. The total acquisition time reported is 40 s for 200 frames.

Fig. 5

B-raster mode acquisition protocol using a SS-OCT system

Figure 6(a) shows a selected B-mode motion frame and Fig. 6(d) shows a 2-D profile of that frame. It is noted that the motion frame is composed of a high spatial frequency varying wave modulated by a low spatial frequency wave. The first wave is produced by the aliasing effect of sampling the propagating shear wave since its velocity is faster than the acquisition speed. The second wave represents the CrW and it is recovered by extracting the envelope of the signal using the Hilbert transform approach [Figs. 6(b) and 6(e)]. Using the phase derivative method described in Sec. 2.3.1, the shear wave speed can be estimated. Figure 6(c) shows the space–time 2-D map at a depth $z = 0.5 mm$ , where the slope indicates the inverse of the shear wave speed. Subsequently, the spatial phase $\tilde{ϕ} (x)$ is estimated [Fig. 6(f)] by using a localized time-domain fitting approach since it is more robust to noise and distortion. Finally, the 2-D shear speed map is calculated using Eq. (18).

Fig. 6

(a) B-mode motion frame of a CrW. (b) Envelope extraction of (a). (c) Space–time map generated from the motion of the envelope in (b). (d) Profile of the B-mode motion frame in (a). (e) Profile of the envelope in (b). (f) Profile of the wave propagation extracted from (c). The color bar scale is in radians related to the OCT phase-shift due to the particle motion in the phantom.

3.3.2.

Swept crawling wave method

This method is analogous to the CrW method; however, the relation $Δ ω ≪ ω$ is no longer required. In this case, the CrW speed is much higher and detecting the slow wave without aliasing is not possible. However, if the sampling theory is applied to Eq. (14), a mathematical model can be found in order to recover the shear wave velocity. The acquisition protocol is based on the formation of M-mode scans for each lateral position; that is, for a given location $x_{0}$ , $M$ depth A-lines are acquired at a rate of 20 kHz (time resolution $Δ_{M} = 50 μ s$ ) as shown in Fig. 7.

Fig. 7

M-mode acquisition protocol for SCrW and SW methods using PhS-OCT. In the particular case in which the acquisition is synchronized with the stimulus, the protocol is called M–B mode used by the SWP method.

Subsequently, the galvo-scanner redirects the laser beam to the next lateral location $x_{1}$ . The M-scanning is performed consecutively for each position $x_{n}$ until the entire field of view of 30 mm is scanned. It should be noted that while the scanning process takes place, the CrW is propagating simultaneously.

The M-mode approach can be modeled by discretizing the time and lateral position such that $t = n Δ_{t}$ , $x = n Δ_{x}$ , for $n = 0, 1, 2, \dots, N - 1$ , and $N$ is the number of A-lines acquired along the 30 mm of lateral field. $Δ_{x}$ is the lateral sampling resolution. For $N = 1000$ , then $Δ_{x} = 30 μ m$ . $Δ_{t}$ is the lag time at each position $x_{n}$ calculated as the sum of the time taken by the system to acquire $M = 300$ A-lines (15 ms) plus the time given to the galvo-scanner (5 ms) to redirect the laser beam to adjacent lateral positions. Then, $Δ_{t} = 20 ms$ . Equation (14) can be rewritten as a function of $Δ_{x}$ and $Δ_{t}$ as

Eq. (22)

{| U (n) |}^{2} = 2 A^{2} e^{- α D} {\cosh (2 α n Δ_{x}) + \cos [Δ ω n Δ_{t} + 2 (k + \frac{Δ k}{2}) n Δ_{x} - φ_{0}]} .

If the ROI is close to the origin, the hyperbolic cosine factor varies slowly and can be considered as a DC factor which can be diminished by using the appropriate filter, leading to

Eq. (23)

{| U (n) |}^{2} = B (\cos {[Δ ω Δ_{t} + 2 (k + \frac{Δ k}{2}) Δ_{x}] n - φ_{0}}),

where

B = 2 A^{2} e^{- α D}

. Then the spatial discrete wave number of the signal is

Δ_{x} k_{SCrW} = Δ ω Δ_{t} + 2 (k + \frac{Δ k}{2}) Δ_{x}

, where

k_{SCrW}

is the spatial continuous wave number, and

Eq. (24)

k_{SCrW} = Δ ω \frac{Δ_{t}}{Δ_{x}} + 2 (k + \frac{Δ k}{2}) .

Equation (24) enables the recovery of

k

(needed for estimating the local shear wave speed) based on

k_{SCrW}

. Moreover, if wave number varies along the lateral axis [

k (x)

], then

k_{SCrW} (x)

can be found using two estimators: local fitting and Hoyt’s algorithm.³⁰ Figure 8(a) shows the 2-D version of the harmonic signal of Eq. (23). It can be noted that

k_{SCrW} (x) > k (x)

; therefore, there are more harmonic cycles along the lateral axis, which may improve the calculation of

k (x)

by using the Hoyt’s estimator (described in Sec. 2.3.2). Figure 8(b) shows a 1-D profile of the 2-D map as well as its filtered version. Finally, the 2-D SWSM can be computed using local fitting and the Hoyt’s estimator.

Fig. 8

(a) Amplitude map generated using the SCrW method. (b) Profile section extracted from the red square region in (a) showing the unfiltered and filtered versions. (c) Amplitude map generated using the SW method. (d) Profile extracted from the red square region in (c). The color bar scale is in radians related to the OCT phase-shift due to the particle motion in the phantom.

3.3.3.

Standing wave method

As described in Sec. 2.2.2, an SW is generated when both harmonic vibrators have the same frequency using the setup configuration in Fig. 3. The frequency of operation used here was $f = 400 Hz$ . We used the same M-mode acquisition protocol (Fig. 7) and parameters as in the swept crawling wave (SCrW) case. Then, Eq. (22) becomes

Eq. (25)

{| U (n) |}^{2} = 2 A^{2} e^{- α D} [\cosh (2 α n Δ_{x}) + \cos (2 k n Δ_{x} - φ_{0})] .

Using the appropriate filter, the DC bias can be suppressed, so that

Eq. (26)

{| U (n) |}^{2} = B [\cos ([2 k Δ_{x}] n - φ_{0})] .

Then the SW number

k_{SW} = 2 k

can be recovered by applying Hoyt’s method³⁰ to the previous signal. However, it is important to note that the number of cycles of

{| U (n) |}^{2}

in Eq. (26) will be much less than that in the SCrW case, which will decrease the effectiveness of the estimation. Figure 8(c) shows the 2-D version of the harmonic signal of Eq. (26), and a 1-D profile version is shown in Fig. 8(d). Finally, the 2-D SWSM is calculated using Eqs. (18) and (20).

3.3.4.

Shear wave propagation method

This method makes use of one vibration source in the setup shown in Fig. 3. A continuous harmonic wave $W_{left} (x, t)$ is sent through the right direction of the phantom at a frequency of $f = 400 Hz$ . The M-mode acquisition is used with the same parameters as in the SCrW and SW methods. The main idea is to generate space–time representations of the propagating shear wave for all depths within the ROI as performed analogously with slow CrWs in the CrW method. The main difficulty lies in the high speed of the shear wave compared to the acquisition speed. If M-mode scanning is performed in the lateral position $x_{0}$ , when the beam moves to the subsequent position $x_{1}$ , the shear wave will propagate several cycles adding an unknown phase to Eq. (7). The sampling time and lateral position can be discretized as $t = m Δ_{M}$ , and $x = n Δ_{x}$ , for $m = 0, 1, 2, \dots, M - 1$ , and $n = 0, 1, 2, \dots, N - 1$ . $Δ_{x}$ and $Δ_{M}$ are the lateral and time sampling resolution, respectively. However, when the beam is redirected from $x_{1}$ to $x_{2}$ , the time required by the galvo-scanner $Δ_{g} = 5 ms$ needs to be added to $t$ . Then $t^{'} = m Δ_{M} + n Δ_{g}$ , and Eq. (7) becomes

Eq. (27)

W (m, n) = A e^{- α (n Δ_{x} + \frac{D}{2})} e^{- j [k (n Δ_{x} + \frac{D}{2}) - ω m Δ_{M}]} e^{j (ω n Δ_{g})},

where the last exponential term adds an extra phase factor

β (n) = ω n Δ_{g}

which is dependent on

n

and will distort the correct generation of space–time representations of the wave propagation. However,

β (n)

can be set to a constant for all values of

n

with the following synchronization requirement:

Eq. (28)

β (n) = β (n + 1) .

Then,

ω Δ_{g} = 2 π l

is needed and

l \in N^{+}

. This restriction is not difficult to satisfy. In fact, using the parameters of the M-mode configuration described in the SCrW method,

ω Δ_{g} = (2 π) (400) (5 \times 10^{- 3}) = 4 π

, and synchronization is ensured. In such a case, the acquisition protocol is called M–B mode as described in Fig. 7. Figure 9(a) shows the space–time map calculated using Eq. (27). The change of slope (related to the shear wave velocity) from one medium to the other is evident. Applying the phase derivative method, the spatial phase

ϕ (x_{n}) = k n Δ_{x}

can be retrieved from Eq. (27) [Fig. 9(b)], and the 2-D shear wave speed may be computed.

Fig. 9

(a) Space–time map describing SWP taken at a depth of 0.8 mm in the phantom. The color bar scale is in radians related to the OCT phase-shift due to the particle motion in the phantom. (b) Profile describing the wave propagation extracted from (a).

3.3.5.

Tone-burst propagation method

Tone-burst propagation (TBP) is the final method compared in this study. It has been widely applied in recent years of OCE research. This method uses one harmonic source to create a tone-burst of an SAW (Rayleigh wave) propagating within the ROI of the phantom. The tone is formed by one cycle of a harmonic wave of $f = 400 Hz$ . The M–B mode approach is used for generating space–time representations of the propagating pulse. In this case, the OCT acquisition and the function generator that produces the tone-burst are triggered at the same time by the computer; thus the synchronization is ensured (Fig. 10). For each lateral position ( $N = 1000$ locations), a tone is burst and $M = 300$ A-lines are acquired. Then the system is able to detect the tone shape in a single M-mode scanning, which requires 15 ms. The time given to the galvo-scanner for redirecting the beam between consecutive lateral positions is $Δ_{g} = 5 ms$ .

Fig. 10

M–B mode acquisition protocol for the TBP method using a PhS-OCT system.

One thousand M-mode maps are acquired using the described protocol in order to cover the ROI lateral axis of 30 mm. A 2-D M-mode map is formed using $300 (time axis) \times 1005 (depth axis)$ elements. We rearranged the scanning into a three-dimensional (3-D) volume of $1000 \times 1005 \times 300$ elements in the lateral, axial, and time axes, respectively. 2-D spatial motion frames at a frame rate of 20 kHz (time resolution $Δ_{M} = 50 μ s$ ) can be obtained from the 3-D volume. Figures 11(a)–11(c) shows the propagation of the tone-burst in three different frames. Moreover, a 2-D profile cut of the 3-D volume at a given depth $z$ will provide a 2-D space–time map [Fig. 11(d)]. Here, it is evident that the slope of the tone trajectory changes from one medium to the other. By applying peak tracking for the time-of flight estimator described in Sec. 2.3.3 on the main peak of the tone in the space–time map [Fig. 11(e)], the Rayleigh wave speed can be estimated. Finally, a 2-D shear wave velocity map is calculated using Eq. (11).

Fig. 11

(a)–(c) B-mode motion frames at different instants depicting the propagation of a tone-burst. (d) Space–time map extracted from the wave propagation using the TBP method taken at a depth of 0.8 mm in the phantom. The color bar scale is in radians related to the OCT phase-shift due to the particle motion in the phantom. (e) Propagation profile extracted from (d).

3.4.

Numerical Simulations

A numerical simulation of elastic wave propagation using the $k$ -space method has been conducted using the k-Wave toolbox in MATLAB (The Mathworks Inc., Natick, Massachusetts) for comparison with experimental results. The numerical grid is depicted in Fig. 12 where the vibration sources, elastic media, boundary conditions, and ROI are indicated. A two-sided media was defined with corresponding shear wave speeds of $c_{s}^{1} = 3.11 m / s$ , and $c_{s}^{2} = 4.78 m / s$ , respectively. The size of the 2-D grid is $78 mm (600 elements) \times 43 mm (300 elements)$ in the lateral and axial directions, respectively. The two vibration sources introduce axial up-and-down particle displacement in the medium in order to generate approximately planar shear waves propagating in the lateral direction. The effect of surface waves was not taken into account in this simulation. Figures 13(a)–13(c) shows three 2-D snapshots of axial direction particle displacement illustrating a TBP pattern at different times within the ROI in order to emulate the same conditions as the PhS-OCT system. Subsequently, we conducted five numerical simulations using one and two vibration sources according to the setup and acquisition protocol described in each dynamic OCE method. For the TBP method, Fig. 13(d) shows the space–time map for a given depth. Here, the slope difference between the two media is evident.

Fig. 12

Numerical simulation setup. Vibration sources, ROI, and medium properties in a 2-D finite grid representing the elastic two-sided medium using the numerical simulator k-Wave toolbox in MATLAB.

Fig. 13

(a)–(c) 2-D snapshots of axial direction particle displacement illustrating a TBP pattern at different times within the ROI using numerical simulations. (d) Space–time map retrieved from the simulated wave motion taken at a depth of 0.8 mm in spatial grid. The color bar scale is in $μ m$ related to the particle displacement.

3.5.

Metrics

Spatial lateral resolution of any OCE method provides important information about the effectiveness of representing the mechanical properties of small elements in a phantom. Assuming that the transition from the 10% to the 15% region in the phantom is a step function, we estimate the resolution by means of a sigmoid function fitted to the estimated shear wave speed profile as described by Rouze et al.³⁵ Then the speed profile is fitted to the function

Eq. (29)

c (x) = c_{10 %} + (c_{15 %} - c_{10 %}) \frac{1}{1 + e^{\frac{x_{0} - x}{τ}}},

where

c_{10 %}

and

c_{15 %}

are the wave velocities in the 10% and 15% regions, respectively,

x_{0}

is the lateral position where the speed transition happens, and

τ

represents the width of the transition. The estimation of

τ

allows the calculation of the resolution

R_{2080}

, defined as the distance from 20% to the 80% of the transition.³⁵ Then the lateral resolution is equivalent to

Eq. (30)

R_{2080} = 2 τ \ln (4) .

The contrast to noise ratio (CNR) parameter provides information of contrast between the two shear wave speed regions (the 10% and 15% concentration regions of the phantom) in the presence of noise related to the OCT acquisition system and the dynamic OCE method. The CNR is given by

Eq. (31)

CNR = \frac{| μ_{15 %} - μ_{10 %} |}{\sqrt{σ_{15 %}^{2} + σ_{10 %}^{2}}},

where

μ_{10 %}

,

μ_{15 %}

are the mean values, and

σ_{10 %}

,

σ_{15 %}

are the standard deviation values of the shear wave speed calculated within a selected rectangular region in the 10% and 15% regions, respectively. The size of these regions were

0.5 \times 4 {mm}^{2}

for the experimental results and

8 \times 8 {mm}^{2}

for the simulated results.

The accuracy error is calculated in terms of how close the average value of speed ( $μ_{10 %}$ or $μ_{15 %}$ ) is in the selected region to the ground truth value ( $μ_{gt}$ ). The percentage accuracy error is given by

Eq. (32)

E_{a} = \frac{| μ_{10 / 15 %} - μ_{gt} |}{μ_{gt}} 100 % .

Finally, the precision is related to the standard deviation of the measured speed in the selected region (

σ_{10 %}

or

σ_{15 %}

). The percentage precision error is calculated as

Eq. (33)

E_{p} = \frac{σ_{10 / 15 %}}{μ_{10 / 15 %}} 100 % .

4. Results

4.1.

Mechanical Measurements

Stress–relaxation curves obtained with the compression test were fitted using a four-parameter FD-SLS model [Fig. 14(a)] by using a nonlinear constrained optimization approach for both phantom concentrations as shown in Fig. 14(b). The FD-SLS model is given by

Eq. (34)

σ + \frac{η}{E_{1}} D_{t}^{α} σ = E_{0} ϵ + η \frac{E_{0} + E_{1}}{E_{1}} D_{t}^{α} ϵ,

where

σ

is the stress,

ϵ

is the strain,

D_{t}^{α}

is the fractional derivative operator, and the four parameters of the model:

E_{0}

(left spring),

E_{1}

(right spring), and

η

,

α

(spring-pot) as described in Schmidt et al.³¹ Once the four parameters are found using the optimization approach, Eq. (34) can be transformed in the frequency domain to

Eq. (35)

\tilde{σ} = \frac{E_{0} + {(i 2 π f)}^{α} η \frac{E_{0} + E_{1}}{E_{1}}}{1 + {(i 2 π f)}^{α} \frac{η}{E_{1}}} \tilde{ϵ} = E^{*} \tilde{ϵ},

where

f

is the frequency,

\tilde{σ}

and

\tilde{ϵ}

are the Fourier domain stress and strain, respectively, and

E^{*}

is the frequency-dependent Young’s modulus. Subsequently, plots of shear wave speed versus excitation frequency [Fig. 14(c)] were obtained by using

E^{*} = 3 G^{*}

in Eq. (9).

Fig. 14

(a) FD-SLS model. (b) Stress relaxation test of the 10% and 15% concentration phantom using the FD-SLS model. (c) Shear wave speed frequency-dependent curves calculated from (b) for both phantom concentrations.

A small dispersion behavior is observed that highlights the low-viscosity component of both phantoms. Then we can neglect the dispersion effect during the wave propagation. We found that the average shear wave velocities of the 10% and 15% gelatin phantoms for a 400-Hz excitation are $3.11 \pm 0.05 m / s$ and $4.78 \pm 0.06 m / s$ , respectively.

4.2.

Simulation Results

A numerical solution of the wave field was created using the k-Wave toolbox in MATLAB using the same dimensions and mechanical properties as the actual phantom being measured. SWSMs that represent the 2-D shear wave speed for each lateral and depth position in a color-coded fashion were obtained by using the five dynamic OCE methods proposed in this study: CrW, SCrW, SW, SWP, and TBP as shown in Fig. 15. For each case, a mean-speed lateral profile was calculated in order to visualize the speed transition from the 10% to the 15% medium. Two noise conditions were simulated: corrupted signals with $SNR = 24 dB$ , and $SNR = 16 dB$ , which corresponded to vibration noise floors of RMS values $0.063 μ m$ and $0.15 μ m$ , respectively. Those noise levels were observed and extracted during the analysis of experimental data in different scenarios. The noise was defined with a zero mean Gaussian distribution and it was added to the particle displacement signal obtained in the numerical simulation. For all cases, we evaluated accuracy, precision of shear wave speed, lateral resolution $R_{2080}$ , and CNR as shown in Table 1. The ground truth speed values were defined in the setup of the simulation described in Sec. 3.4.

Fig. 15

Spatially resolved 2-D SWSM and 1-D profiles obtained from all OCE methods during numerical simulations: (a) CrW using phase derivative estimator, (b) SW using Hoyt’s estimator, (c) SCrW using local fitting estimator, (d) SCrW using Hoyt’s estimator, (e) SWP using phase derivative estimator, and (f) TBP using time-of-flight estimator. Color bar scale is in $m / s$ .

Table 1

Accuracy error and precision of shear wave speed, CNR, and resolution (R2080), for all methods extracted from numerical simulations under different levels of noise corruption. PhD, phase derivative estimator; Hoyt, Hoyt’s estimator; TofF, time of flight approach. Ground truth values: cshear 10%=3.11±0.05 m/s and cshear 15%=4.78±0.06 m/s.

		Gel 10%		Gel 15%		CNR	R2080 (mm)
		Accuracy error (%)	Precision error (%)	Accuracy error (%)	Precision error (%)	CNR	R2080 (mm)
CrW	PhD—SNR 24 dB	0.69	6.23	11.67	14.39	2.78	0.73
CrW	PhD—SNR 16 dB	0.74	6.27	11.73	14.51	2.76	0.81
\|SCrW	Fitting—SNR 24 dB	0.22	7.43	13.79	18.18	2.29	0.42
	Fitting—SNR 16 dB	2.93	13.55	21.66	21.32	1.99	3.38
	Hoyt—SNR 24 dB	11.22	11.03	14.38	19.47	1.78	0.55
	Hoyt—SNR 16 dB	55.10	19.97	74.64	27.27	1.43	4.29
SW	Hoyt—SNR 24 dB	1.87	1.10	31.54	3.36	14.57	1.33
SW	Hoyt—SNR 16 dB	1.88	1.10	31.80	4.34	11.36	1.34
SWP	PhD—SNR 24 dB	0.08	2.47	9.42	2.10	15.80	0.69
SWP	PhD—SNR 16 dB	3.55	14.60	15.80	13.75	2.59	0.70
TBP	TofF—SNR 24 dB	0.87	3.74	1.16	1.31	12.74	1.20
TBP	TofF—SNR 15 dB	1.25	3.83	1.42	1.68	11.67	1.25

4.3.

Experimental Results

SWSMs were obtained by using the proposed dynamic OCE methods (CrW, SCrW, SW, SWP, and TBP) as shown in Fig. 16. For each case, a mean-speed lateral profile was calculated in order to visualize the speed transition from the 10% to the 15% medium. The accuracy and precision errors were calculated within a square region in each side of the SWSM representing the 10% and 15% gelatin regions of the phantom. Here, the mechanical measurement results were used as the ground truth for all methods. Mean speed and standard deviation for both regions using all techniques can be found in Table 2, whereas accuracy and precision error values are shown in Table 3. Figure 17(a) shows a bar graph of the estimated shear wave speeds for all methods including mechanical measurements. Figure 17(b) shows a normalized sigmoid function after being fitted to the actual speed profiles for all techniques. $R_{2080}$ is also calculated and shown in Table 3. The axial resolution of SWSM is the same for all techniques and is equal to the optical axial resolution of the OCT system ( $10 μ m$ ), which currently limits performance in this metric. CNR values were also calculated for all methods and are shown in Table 3. Finally, the last metric used in this comparative study is acquisition time. Time is important for in vivo applications since the target subjected to measurement should be static during the acquisition process. For all techniques, the acquisition time was measured and is reported in Table 3.

Fig. 16

Spatially resolved 2-D SWSM and 1-D profiles obtained from all OCE methods during experiments: (a) CrW using phase derivative estimator, (b) SW using Hoyt’s estimator, (c) SCrW using local fitting estimator, (d) SCrW using Hoyt’s estimator, (e) SWP using phase derivative estimator, and (f) TBP using time-of-flight estimator. Color bar scale is in $m / s$ .

Table 2

Mean and standard deviation of shear wave velocity estimation for all methods in both phantom gelatin concentrations. Ground truth values: cshear 10%=3.11±0.05 m/s and cshear 15%=4.78±0.06 m/s.

	Gel 10%		Gel 15%
	Mean (m/s)	Standard deviation (m/s)	Mean (m/s)	Standard deviation (m/s)
CrW	3.49	0.15	4.86	0.19
SCrW	3.87	0.37	5.36	0.62
SW	3.87	0.52	6.57	1.06
SWP	3.30	0.30	5.01	0.14
TBP	3.11	0.21	4.59	0.33

Table 3

Accuracy error and precision of shear wave speed estimation, CNR, resolution (R2080), and acquisition time for all methods extracted from experimental results. Ground truth values: cshear 10%=3.11±0.05 m/s and cshear 15%=4.78±0.06 m/s.

	Gel 10%		Gel 15%		CNR	R2080 (mm)	Time (s)
	Accuracy error (%)	Precision error (%)	Accuracy error (%)	Precision error (%)	CNR	R2080 (mm)	Time (s)
CrW	12.22	3.14	1.67	3.97	5.66	0.99	40
SCrW	24.44	7.74	12.13	12.97	2.06	1.36	20
SW	24.44	10.88	37.45	22.18	2.29	1.30	20
SWP	6.11	6.28	4.81	2.93	5.17	0.82	20
TBP	0.16	4.39	3.97	6.90	3.80	0.88	20

Fig. 17

(a) Comparison of average shear wave speed for all methods in both phantom concentrations. (b) Normalized sigmoid curves describing the transition from 10% to 15% medium in all methods. Ground truth values: $c_{shear 10 %} = 3.11 \pm 0.05 m / s$ , and $c_{shear 15 %} = 4.78 \pm 0.06 m / s$ .

An SWSM of a chicken tibialis anterior muscle with a thermal lesion was obtained by using the TBP method (Fig. 18) since it reported the highest accuracy in the study with phantoms. The average shear wave speeds within a rectangular region located inside and outside the lesion are $5.71 \pm 0.82 m / s$ and $7.47 \pm 0.52 m / s$ , respectively. The CNR was calculated to be 1.85.

Fig. 18

Elastography of a chicken tibialis anterior muscle. (a) B-mode image depicting the thermal lesion (gamma correction value 0.3). (b) SWSM using TBP method. Color bar scale is in $m / s$ .

5. Discussion and Conclusion

In this study, five dynamic OCE methods were employed to estimate the shear wave speed of a two-sided gelatin tissue-mimicking phantom (10% and 15% concentration): CrW, SCrW, SW, SWP, and TBP. As shown in Table 1, numerical results predict that the TBP method has the most accurate speed estimation on both sides of the phantom ( $accuracy error < 1.5 %$ ) for different noise conditions. This outcome is then confirmed experimentally as reported in Table 3 for the TBP method ( $accuracy error < 4 %$ ). While SWP has the second most accurate results ( $accuracy error < 6.2 %$ ) in experiments, it also provides the second best CNR. The remaining methods overestimate the shear wave speed on both sides of the phantom in simulations and experiments. In particular, the SW experimental result is the least accurate with the most variable estimation since the nature of this approach relies on the recovery of $k (x)$ using just six cycles of the signal ${| U (n) |}^{2}$ . SCrW improves this estimation by recovering $k (x)$ using 85 cycles of the same signal confirming what was predicted in simulation. The CrW method reports the best CNR among all cases in simulation and experiments.

There is a marked difference between OCE methods that use one versus two wave vibration sources. TBP and SWP remain more accurate than CrW, SCrW, and SW. This behavior is attributed to the multiple wave reflections at all boundaries of the phantom that are more prevalent in the dual source and steady-state techniques. OCE methods using two sources have twice the number of reflections than one-source methods in experiments. While interference due to reflections is present on both sides of the phantom for the two-source methods [see Figs. 16(a)–16(d)], SWP depicts interference just on the 10% side [Fig. 16(e)]. In addition, TBP is less sensitive to reflections than all other methods because the energy of the propagating pulse decays rapidly compared to a continuous wave. However, a shear wave traveling from one medium to another of different elastic modulus produces a reflection than cannot be avoided. Figure 15(f) shows this interference in the transition between the 10% and 15% media using the TB method during simulations. Experimental results confirm this effect in the transition between the two media for the TBP method [see Fig. 16(f)]. Directional filters can be used in order to mitigate the reflections, especially in samples with prominent boundaries and high gradients of elasticity as proposed by Song et al.¹⁷

Lateral resolution is found to be similar ( $R_{2080} < 1 mm$ ) for CrW, SWP, and TBP; as well as ( $R_{2080} < 1.4 mm$ ) for SCrW and SW during experiments. Figure 17(b) shows the sigmoid curves for all cases depicting the transition of velocity from the 10% to the 15% side. This measurement is sensitive to the filter size (phase derivative method), and window size (Hoyt’s method). There is a compromise between lateral resolution and the estimation technique for each case. According to numerical simulations, the phase derivative estimator reports a better lateral resolution compared to the Hoyt’s estimator, whereas the variability of the former is higher than the latter as shown in Figs. 15(a) and 15(b). Similarly, the time of flight approach used in the TBP method has a higher lateral resolution than the Hoyt’s estimator. This is consistent with experimental results in which SCrW and SW, which utilize the Hoyt’s estimator, have the lowest lateral resolution among all methods. In particular, the CrW method reports the highest lateral resolution during simulation. In experiments, although CrW is grouped into the highest lateral resolution cluster, SWP has the best lateral resolution compared to all cases. The discrepancy in lateral resolution between experimental and simulated results can be attributed to some extent to the mixing process between the 10% and 15% phantom regions. While the phantoms were created in a two-step process to minimize a mixing/bleeding behavior and the imaging was conducted within a half day of manufacturing the phantom, the migration of fluid from the higher concentration region to the lower one was inevitable over time. Therefore, the transition from the 10% to 15% regions was not an ideal step function and may have led to the underestimation of the lateral resolution.

Resolution results are comparable to those obtained using USE. According to the work of Ormachea et al.,³⁶ lateral resolutions in the approximate range of 2 to 4 mm are reported for excitation frequencies up to 500 Hz using the CrW method. Resolution of wave-based elastography techniques is related primarily to the excitation wavelength, and, secondarily, to the resolution of the imaging modality. In our study, the excitation frequency used is 400 Hz, which is comparable to the frequency range used in the US study.³⁶ This explains why lateral resolution in both imaging modalities is in the same order of magnitude. However, OCT is sensitive enough to detect waves with smaller wavelengths produced by excitations of higher frequencies (in the order of kilohertz), which may represent an advantage in resolution compared to US modalities.

Experimental results in a chicken tibialis muscle using the TBP method report shear wave speed values consistent with a previous elastography study in chicken using US.³⁷ Here, the estimated shear modulus of the left side of a chicken tibialis anterior muscle of $25.3 \pm 2.2 kPa$ can be translated to a shear wave speed of $5.03 \pm 0.22 m / s$ using Eq. (6), which is similar to our results for normal muscle ( $5.71 \pm 0.82 m / s$ ). The thermal lesion is stiffer when compared to normal muscle tissue, which enables its detection as shown in Fig. 8.

The acquisition time is an important factor to be considered when using an OCE method for in vivo applications. CrW require twice the time of the other methods during the acquisition process. In CrW, time can be reduced by scanning fewer frames at the cost of losing accuracy during the estimation procedure. In TBP, the acquisition time can be reduced to 15 s at the same lateral sampling resolution if we ensure that the burst peak can be captured by using $M = 200$ time sampling repetitions. For all methods, the time can be reduced by reducing the number of samples taken in the lateral direction. In addition, the acquisition speed can be dramatically increased by using an OCT system with a faster acquisition A-line rate as proposed in Singh et al.,³⁸ and Song et al.³⁹ In addition, holographic techniques⁴⁰^,⁴¹ demonstrated the fast acquisition of surface wave propagation in tissue, which offers an alternative to implement those OCE methods in real time.

In summary, all five evaluated dynamic OCE methods have both advantages and limitations. In particular, the TBP method reports the best accuracy on both sides of the phantom during experiments and numerical simulations. It is also less sensitive to reflections which makes it a promising technique for the quantitative estimation of elasticity in biotissue. The dispersion effect of the wave speed versus excitation frequency is an important behavior caused by viscoelastic media: however, this study only considered elastic media. Future work will include viscoelastic materials and we will investigate the frequency-dependent estimation of shear wave speed in order to calculate the viscosity parameter, which has been demonstrated to be an important biomarker or the characterization of tissue such as human liver⁴² and human cornea.⁴³ In addition, material for further study includes: (a) special cases such as depth-layered media and hard inclusions in softer background, and (b) how loading conditions using other technologies can be used for the better characterization of elasticity. In particular, the work of Mohan and Oldenburg⁴⁴ offers the possibility to resolve depth-layered media using Rayleigh wave propagation by solving a constrained inverse problem.

Disclosures

No conflicts of interest, financial or otherwise, are declared by the authors.

Acknowledgments

The instrumentation engineering development of this research benefitted from support of the II-VI Foundation. Fernando Zvietcovich is supported by the Fulbright Program (U.S.A. Department of State), and Fondo para la Innovacion, la Ciencia, y la Tecnologia (FINCyT) 097-FINCYT-BDE-2014 (Peruvian Government).

References

1.

L. Sandrin et al., “Transient elastography: a new noninvasive method for assessment of hepatic fibrosis,” Ultrasound Med. Biol., 29 (12), 1705 –1713 (2003). http://dx.doi.org/10.1016/j.ultrasmedbio.2003.07.001 USMBA3 0301-5629 Google Scholar

2.

C. Schmitt et al., “Noninvasive vascular elastography: toward a complementary characterization tool of atherosclerosis in carotid arteries,” Ultrasound Med. Biol., 33 (12), 1841 –1858 (2007). http://dx.doi.org/10.1016/j.ultrasmedbio.2007.05.020 USMBA3 0301-5629 Google Scholar

3.

J. Rogowska et al., “Optical coherence tomographic elastography technique for measuring deformation and strain of atherosclerotic tissues,” Heart, 90 (5), 556 –562 (2004). http://dx.doi.org/10.1136/hrt.2003.016956 Google Scholar

4.

B. P. Donald et al., “Visualization and quantification of breast cancer biomechanical properties with magnetic resonance elastography,” Phys. Med. Biol., 45 (6), 1591 –1610 (2000). http://dx.doi.org/10.1088/0031-9155/45/6/314 PHMBA7 0031-9155 Google Scholar

5.

E. T. Detorakis and I. G. Pallikaris, “Ocular rigidity: biomechanical role, in vivo measurements and clinical significance,” Clin. Exp. Ophthalmol., 41 (1), 73 –81 (2013). http://dx.doi.org/10.1111/ceo.2013.41.issue-1 Google Scholar

6.

K. J. Parker, M. M. Doyley and D. J. Rubens, “Imaging the elastic properties of tissue: the 20 year perspective,” Phys. Med. Biol., 56 (1), R1 (2011). http://dx.doi.org/10.1088/0031-9155/56/1/R01 PHMBA7 0031-9155 Google Scholar

7.

R. Muthupillai et al., “Magnetic resonance elastography by direct visualization of propagating acoustic strain waves,” Science, 269 (5232), 1854 –1857 (1995). http://dx.doi.org/10.1126/science.7569924 SCIEAS 0036-8075 Google Scholar

8.

W. Drexler et al., “Optical coherence tomography today: speed, contrast, and multimodality,” J. Biomed. Opt., 19 (7), 071412 (2014). http://dx.doi.org/10.1117/1.JBO.19.7.071412 JBOPFO 1083-3668 Google Scholar

9.

J. A. Mulligan et al., “Emerging approaches for high-resolution imaging of tissue biomechanics with optical coherence elastography,” IEEE J. Sel. Top. Quantum Electron., 22 (3), 246 –265 (2016). http://dx.doi.org/10.1109/JSTQE.2015.2481705 IJSQEN 1077-260X Google Scholar

10.

J. M. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express, 3 (6), 199 –211 (1998). http://dx.doi.org/10.1364/OE.3.000199 OPEXFF 1094-4087 Google Scholar

11.

R. K. Wang, Z. Ma and S. J. Kirkpatrick, “Tissue Doppler optical coherence elastography for real time strain rate and strain mapping of soft tissue,” Appl. Phys. Lett., 89 (14), 144103 (2006). http://dx.doi.org/10.1063/1.2357854 APPLAB 0003-6951 Google Scholar

12.

C. Li et al., “Determining elastic properties of skin by measuring surface waves from an impulse mechanical stimulus using phase-sensitive optical coherence tomography,” J. R. Soc. Interface, 9 (70), 831 –841 (2012). http://dx.doi.org/10.1098/rsif.2011.0583 1742-5689 Google Scholar

13.

R. Manapuram et al., “Estimation of shear wave velocity in gelatin phantoms utilizing PhS-SSOCT,” Laser Phys., 22 (9), 1439 –1444 (2012). http://dx.doi.org/10.1134/S1054660X12090101 LAPHEJ 1054-660X Google Scholar

14.

S. Song, Z. Huang and R. K. Wang, “Tracking mechanical wave propagation within tissue using phase-sensitive optical coherence tomography: motion artifact and its compensation,” J. Biomed. Opt., 18 (12), 121505 (2013). http://dx.doi.org/10.1117/1.JBO.18.12.121505 JBOPFO 1083-3668 Google Scholar

15.

C. Li et al., “Full skin quantitative optical coherence elastography achieved by combining vibration and surface acoustic wave methods,” Proc. SPIE, 9322 93220O (2015). http://dx.doi.org/10.1117/12.2075666 PSISDG 0277-786X Google Scholar

16.

A. Nahas et al., “From supersonic shear wave imaging to full-field optical coherence shear wave elastography,” J. Biomed. Opt., 18 (12), 121514 (2013). http://dx.doi.org/10.1117/1.JBO.18.12.121514 JBOPFO 1083-3668 Google Scholar

17.

S. Song et al., “Quantitative shear-wave optical coherence elastography with a programmable phased array ultrasound as the wave source,” Opt. Lett., 40 (21), 5007 –5010 (2015). http://dx.doi.org/10.1364/OL.40.005007 OPLEDP 0146-9592 Google Scholar

18.

S. Wang and K. V. Larin, “Noncontact depth-resolved micro-scale optical coherence elastography of the cornea,” Biomed. Opt. Express, 5 (11), 3807 –3821 (2014). http://dx.doi.org/10.1364/BOE.5.003807 BOEICL 2156-7085 Google Scholar

19.

J. Li et al., “Dynamic optical coherence tomography measurements of elastic wave propagation in tissue-mimicking phantoms and mouse cornea in vivo,” J. Biomed. Opt., 18 (12), 121503 (2013). http://dx.doi.org/10.1117/1.JBO.18.12.121503 JBOPFO 1083-3668 Google Scholar

20.

S. Wang and K. V. Larin, “Shear wave imaging optical coherence tomography (SWI-OCT) for ocular tissue biomechanics,” Opt. Lett., 39 (1), 41 –44 (2014). http://dx.doi.org/10.1364/OL.39.000041 OPLEDP 0146-9592 Google Scholar

21.

A. Ahmad et al., “Magnetomotive optical coherence elastography using magnetic particles to induce mechanical waves,” Biomed. Opt. Express, 5 (7), 2349 –2361 (2014). http://dx.doi.org/10.1364/BOE.5.002349 BOEICL 2156-7085 Google Scholar

22.

G. F. Pinton, J. J. Dahl and G. E. Trahey, “Rapid tracking of small displacements with ultrasound,” IEEE Trans. Ultrasonics Ferroelectrics Freq. Control, 53 (6), 1103 –1117 (2006). http://dx.doi.org/10.1109/TUFFC.2006.1642509 Google Scholar

23.

T. Loupas, R. B. Peterson and R. W. Gill, “Experimental evaluation of velocity and power estimation for ultrasound blood flow imaging, by means of a two-dimensional autocorrelation approach,” IEEE Trans. Ultrasonics Ferroelectrics Freq. Control, 42 (4), 689 –699 (1995). http://dx.doi.org/10.1109/58.393111 Google Scholar

24.

K. F. Graff, Wave Motion in Elastic Solids, Dover Publications(2012). Google Scholar

25.

E. C. Elegbe and S. A. McAleavey, “Single tracking location methods suppress speckle noise in shear wave velocity estimation,” Ultrasonic Imaging, 35 (2), 109 –125 (2013). http://dx.doi.org/10.1177/0161734612474159 ULIMD4 0161-7346 Google Scholar

26.

T. Glozman and H. Azhari, “A method for characterization of tissue elastic properties combining ultrasonic computed tomography with elastography,” J. Ultrasound Med., 29 (3), 387 –398 (2010). http://dx.doi.org/10.7863/jum.2010.29.3.387 JUMEDA 0278-4297 Google Scholar

27.

F. Zvietcovich et al., “Experimental classification of surface waves in optical coherence elastography,” Proc. SPIE, 9710 97100Z (2016). http://dx.doi.org/10.1117/12.2211420 PSISDG 0277-786X Google Scholar

28.

Z. Wu et al., “Sonoelastographic imaging of interference patterns for estimation of the shear velocity of homogeneous biomaterials,” Phys. Med. Biol., 49 (6), 911 –922 (2004). http://dx.doi.org/10.1088/0031-9155/49/6/003 PHMBA7 0031-9155 Google Scholar

29.

Z. Hah et al., “Integration of crawling waves in an ultrasound imaging system. Part 2: signal processing and applications,” Ultrasound Med. Biol., 38 (2), 312 –323 (2012). http://dx.doi.org/10.1016/j.ultrasmedbio.2011.10.014 USMBA3 0301-5629 Google Scholar

30.

K. Hoyt, B. Castaneda and K. J. Parker, “Two-dimensional sonoelastographic shear velocity imaging,” Ultrasound Med. Biol., 34 (2), 276 –288 (2008). http://dx.doi.org/10.1016/j.ultrasmedbio.2007.07.011 USMBA3 0301-5629 Google Scholar

31.

A. Schmidt and L. Gaul, “Experimental investigation and numerical treatment of viscoelastic materials,” Proc. Int. Model Anal. Conf., 3 (1), 1557 –1566 (2008). Google Scholar

32.

J. Yao et al., “Angular scan optical coherence tomography imaging and metrology of spherical gradient refractive index preforms,” Opt. Express, 23 (5), 6428 –6443 (2015). http://dx.doi.org/10.1364/OE.23.006428 OPEXFF 1094-4087 Google Scholar

33.

P. Meemon, K. Lee and J. Rolland, “Doppler imaging with dual-detection full-range frequency domain optical coherence tomography,” Biomed. Opt. Express, 1 (2), 537 –552 (2010). http://dx.doi.org/10.1364/BOE.1.000537 BOEICL 2156-7085 Google Scholar

34.

P. Meemon et al., “Crawling wave optical coherence elastography,” Opt. Lett., 41 (5), 847 –850 (2016). http://dx.doi.org/10.1364/OL.41.000847 OPLEDP 0146-9592 Google Scholar

35.

N. C. Rouze et al., “Parameters affecting the resolution and accuracy of 2-D quantitative shear wave images,” IEEE Trans. Ultrasonics Ferroelectrics Freq. Control, 59 (8), 1729 –1740 (2012). http://dx.doi.org/10.1109/TUFFC.2012.2377 Google Scholar

36.

J. Ormachea et al., “Shear wave speed measurements using crawling wave sonoelastography and single tracking location shear wave elasticity imaging for tissue characterization,” IEEE Trans. Ferroelectrics Freq. Control, 63 (9), 1351 –1360 (2016). http://dx.doi.org/10.1109/TUFFC.2016.2576962 Google Scholar

37.

T. K. Koo et al., “Relationship between shear elastic modulus and passive muscle force: an ex-vivo study,” J. Biomechan., 46 (12), 2053 –2059 (2013). http://dx.doi.org/10.1016/j.jbiomech.2013.05.016 Google Scholar

38.

M. Singh et al., “Phase-sensitive optical coherence elastography at 1.5 million A-Lines per second,” Opt. Lett., 40 (11), 2588 –2591 (2015). http://dx.doi.org/10.1364/OL.40.002588 OPLEDP 0146-9592 Google Scholar

39.

S. Song et al., “Strategies to improve phase-stability of ultrafast swept source optical coherence tomography for single shot imaging of transient mechanical waves at 16 kHz frame rate,” Appl. Phys. Lett., 108 (19), 191104 (2016). http://dx.doi.org/10.1063/1.4949469 APPLAB 0003-6951 Google Scholar

40.

C. H. Liu et al., “Non-contact single shot elastography using line field low coherence holography,” Biomed. Opt. Express, 7 (8), 3021 –3031 (2016). http://dx.doi.org/10.1364/BOE.7.003021 BOEICL 2156-7085 Google Scholar

41.

S. Li et al., “Toward soft-tissue elastography using digital holography to monitor surface acoustic waves,” J. Biomed. Opt., 16 (11), 116005 (2011). http://dx.doi.org/10.1117/1.3646211 JBOPFO 1083-3668 Google Scholar

42.

K. R. Nightingale et al., “Derivation and analysis of viscoelastic properties in human liver: impact of frequency on fibrosis and steatosis staging,” IEEE Trans. Ultrasonics Ferroelectrics Freq. Control, 62 (1), 165 –175 (2015). http://dx.doi.org/10.1109/TUFFC.2014.006653 Google Scholar

43.

G. Lombardo et al., “Analysis of the viscoelastic properties of the human cornea using Scheimpflug imaging in inflation experiment of eye globes,” PLoS One, 9 (11), e112169 (2014). http://dx.doi.org/10.1371/journal.pone.0112169 POLNCL 1932-6203 Google Scholar

44.

K. Mohan and A. Oldenburg, “Elastography of soft materials and tissues by holographic imaging of surface acoustic waves,” Opt. Express, 20 18887 –18897 (2012). http://dx.doi.org/10.1364/OE.20.018887 OPEXFF 1094-4087 Google Scholar

Biographies for the authors are not available.

Citation Download Citation

Fernando Zvietcovich, Jannick P. Rolland, Jianing Yao, Panomsak Meemon, and Kevin J. Parker "Comparative study of shear wave-based elastography techniques in optical coherence tomography," Journal of Biomedical Optics 22(3), 035010 (30 March 2017). https://doi.org/10.1117/1.JBO.22.3.035010

Received: 18 November 2016; Accepted: 15 March 2017; Published: 30 March 2017

Access the abstract

JOURNAL ARTICLE
17 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY