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

**Fernando Zvietcovich**,

**Kevin J. Parker**

University of Rochester, Department of Electrical and Computer Engineering, Rochester, New York, United States

**Jannick P. Rolland**,

**Jianing Yao**

University of Rochester, The Institute of Optics, Rochester, New York, United States

**Panomsak Meemon**

University of Rochester, The Institute of Optics, Rochester, New York, United States

Suranaree University of Technology, School of Physics, Institute of Science, Nakhon Ratchasima, Thailand

*J. Biomed. Opt*. 22(3), 035010 (Mar 30, 2017). doi:10.1117/1.JBO.22.3.035010

#### Open Access

## Abstract

**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.

Elastic properties of tissue provide valuable information regarding pathological processes such as liver fibrosis,^{1} arteriosclerosis,^{2}^{,}^{3} breast cancer,^{4} and glaucoma,^{5} 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.^{6} Magnetic resonance elastography, on the other hand, has a wider scanning range ($\u223c10$ to 100 mm) with resolutions on the order of $\u223c10\u2009\u2009mm$.^{7} 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.^{8} 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.^{9}

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).^{9} Quasi-static OCE was first demonstrated by Schmitt.^{10} It is also important to mention the work of Wang et al.^{11} 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.^{12} 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,^{12}^{–}^{15} acoustic radiation force,^{16}^{,}^{17} air-puff excitation,^{18}^{–}^{20} and magnetic transducers^{21} were utilized in the wave generation. In particular, Manapuram et al.^{13} used a swept-source OCT (SS-OCT) for measuring shear wave pulse speed in gelatin phantoms. Similarly, Song et al.^{14} 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.^{15} used the vibration method and surface acoustic waves (SAW) to determine the biomechanical properties of epidermal layers of skin. Nahas et al.^{16} 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^{18} 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.^{21}

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.

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.

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 $\Delta \varphi (z)=\varphi (z,t1)\u2212\varphi (z,t0)$ at two consecutive instants $t0$ and $t1$ ($t0<t1$) for a given lateral position $x0$ is related to the particle displacement in the axial direction by

PhS-OCT possesses a high displacement sensitivity, on the order of nanometers,^{11} 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 $uz(z)$ created by the motion of the air–medium interface can be compensated as discussed in the work of Song et al.^{14}

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

^{22}we chose the Loupas et al.

^{23}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, $\Delta \varphi (z)$ can be represented as

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.

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 $xz$-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

^{24}Then, at a given time $t0$, the particle displacement generated by the shear wave [Fig. 2(a)] can be measured by a PhS-OCT system using Eq. (1), and $uz(x,z)=W(x,z;t0)$. The phase velocity of the shear wave, also referred to as the speed of the shear wave, is given by

^{6}

In general, biological tissues are nearly incompressible and the Poisson’s ratio can be approximated to $\nu =0.5$.^{6} Then the shear modulus $G$ can be related to the elastic modulus (Young’s modulus) $E$ as $E=3G$.^{6} 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 $cshear$ is reduced to the estimation of $\lambda $ or, equivalently, by estimating the wave number $k=2\pi /\lambda $. In addition, the shear wave speed can also be estimated by using the kinematic relationship

^{25}Shear wave velocity in tissue typically ranges from 0.1 to $10\u2009\u2009m/s$;

^{26}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.^{27} 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 $\nu =0.5$ is approximately given by^{27}

^{27}

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 $Wleft$ and $Wright$ parallel to the $x$-axis, and with different frequencies $\omega $ and $\omega +\Delta \omega $, respectively. For a selected depth $z0$, we can describe the propagation of those waves in the $x-t$ domain as

^{28}assuming no reflections in the medium, the square amplitude of the interference $U(x,t)=Wleft+Wright$ is given by

^{28}[Fig. 2(c)]. The phase term $\phi 0=\Delta k2D+\Delta \phi $ is constant and $\Delta \phi =\phi r\u2212\phi l$. The CrW speed $cCrW$ is related to the shear wave speed $cshear$ by

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 $cshear$ from the displacement field measured by PhS-OCT will be presented and discussed in this section.

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.^{29} For a given depth $z0$, 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 $\Delta \omega $ in order to remove the nonharmonic cosh term. Then, assuming $\Delta \omega \u226a\omega $, it is also true that $\Delta k\u226ak$, and

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

Therefore, the spatial wave number can be estimated as

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)=[xn,xn+1,xn+2,\u2026,xn+w\u22121]$, and $q(n)=[tn,tn+1,tn+2,\u2026,tn+w\u22121]$, respectively; and the laterally resolved group speed is estimated as

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.

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 $\u223c1.35$, which is close to that of real tissue. The shape of the phantom is a rectangular cuboid of dimensions $15\xd77\xd73\u2009\u2009cm3$ [Fig. 3(a)].

The Young’s modulus of each side of the phantom was measured by taking three cylindrical samples ($n=3$) $\u223c45-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\u2009\u2009mm/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.^{31} Finally, using Eq. (9) for a material density of $1\u2009\u2009g/cm3$, 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.

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/e2$ 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\u2009\u2009MSamples/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.^{32} The imaging depth of the system was measured to be 5 mm in air ($\u221210\u2009\u2009dB$ sensitive fall-off). The optical lateral resolution was $\u223c30\u2009\u2009\mu m$, and the FWHM of the axial point spread function after dispersion compensation was $10\u2009\u2009\mu 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 ($\Delta \varphi err$) while imaging a static structure.^{33} Result shows $\Delta \varphi err=4.6\u2009\u2009mrad$ when using the Loupas’ algorithm of Eq. (6). The displacement sensitivity is measured as the minimum detectable axial displacement ($uz,min$) by plugging $\Delta \varphi err$ in Eq. (1). We found a displacement of $uz,min=0.358\u2009\u2009nm$. Finally, the maximum axial displacement supported by the system without unwrapping the phase-shift signal ($\Delta \varphi max=\pi $) is $uz,max=0.24\u2009\u2009\mu 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 $\u223c9\u2009\u2009cm$. 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 $\u223c8\u2009\u2009cm$ 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\u2009\u2009\mu 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 $uz,max=0.24\u2009\u2009\mu m$ without unwrapping the Doppler phase signal. Since waves travel $\u223c6\u2009\u2009cm$ 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.

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)].

**F4 :**

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.

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.^{28} The same idea was applied for the first time in OCE as reported by Meemon et al.^{34} The two harmonic actuators are set at frequencies $fleft=400\u2009\u2009Hz$ and $fright=400.4\u2009\u2009Hz$, respectively, yielding a $\Delta f=0.4\u2009\u2009Hz$ or $\Delta \omega \u22482.51\u2009\u2009rad/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\u2009\u2009frames/s$ using Loupas’ approach^{23} 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\u2009\u2009\mu m$, as described in Fig. 5. The total acquisition time reported is 40 s for 200 frames.

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\u2009\u2009mm$, where the slope indicates the inverse of the shear wave speed. Subsequently, the spatial phase $\varphi \u02dc(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).

**F6 :**

(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.

This method is analogous to the CrW method; however, the relation $\Delta \omega \u226a\omega $ 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 $x0$, $M$ depth A-lines are acquired at a rate of 20 kHz (time resolution $\Delta M=50\u2009\u2009\mu s$) as shown in Fig. 7.

Subsequently, the galvo-scanner redirects the laser beam to the next lateral location $x1$. The M-scanning is performed consecutively for each position $xn$ 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\Delta t$, $x=n\Delta x$, for $n=0,1,2,\u2026,N\u22121$, and $N$ is the number of A-lines acquired along the 30 mm of lateral field. $\Delta x$ is the lateral sampling resolution. For $N=1000$, then $\Delta x=30\u2009\u2009\mu m$. $\Delta t$ is the lag time at each position $xn$ 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, $\Delta t=20\u2009\u2009ms$. Equation (14) can be rewritten as a function of $\Delta x$ and $\Delta t$ as

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

^{30}Figure 8(a) shows the 2-D version of the harmonic signal of Eq. (23). It can be noted that $kSCrW(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.

**F8 :**

(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.

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\u2009\u2009Hz$. We used the same M-mode acquisition protocol (Fig. 7) and parameters as in the swept crawling wave (SCrW) case. Then, Eq. (22) becomes

^{30}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).

This method makes use of one vibration source in the setup shown in Fig. 3. A continuous harmonic wave $Wleft(x,t)$ is sent through the right direction of the phantom at a frequency of $f=400\u2009\u2009Hz$. 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 $x0$, when the beam moves to the subsequent position $x1$, 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\Delta M$, and $x=n\Delta x$, for $m=0,1,2,\u2026,M\u22121$, and $n=0,1,2,\u2026,N\u22121$. $\Delta x$ and $\Delta M$ are the lateral and time sampling resolution, respectively. However, when the beam is redirected from $x1$ to $x2$, the time required by the galvo-scanner $\Delta g=5\u2009\u2009ms$ needs to be added to $t$. Then $t\u2032=m\Delta M+n\Delta g$, and Eq. (7) becomes

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\u2009\u2009Hz$. 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 $\Delta g=5\u2009\u2009ms$.

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\u2009(time axis)\xd71005\u2009(depth axis)$ elements. We rearranged the scanning into a three-dimensional (3-D) volume of $1000\xd71005\xd7300$ elements in the lateral, axial, and time axes, respectively. 2-D spatial motion frames at a frame rate of 20 kHz (time resolution $\Delta M=50\u2009\u2009\mu 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).

**F11 :**

(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).

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 $cs1=3.11\u2009\u2009m/s$, and $cs2=4.78\u2009\u2009m/s$, respectively. The size of the 2-D grid is $78\u2009\u2009mm\u2009(600\u2009elements)\xd743\u2009\u2009mm\u2009(300\u2009elements)$ 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.

**F13 :**

(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 $\mu m$ related to the particle displacement.

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.^{35} Then the speed profile is fitted to the function

^{35}Then the lateral resolution is equivalent to

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

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

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

^{31}Once the four parameters are found using the optimization approach, Eq. (34) can be transformed in the frequency domain to

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\xb10.05\u2009\u2009m/s$ and $4.78\xb10.06\u2009\u2009m/s$, respectively.

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\u2009\u2009dB$, and $SNR=16\u2009\u2009dB$, which corresponded to vibration noise floors of RMS values $0.063\u2009\u2009\mu m$ and $0.15\u2009\u2009\mu 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 $R2080$, 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.

**F15 :**

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$.

Gel 10% | Gel 15% | CNR | $R2080$ (mm) | ||||
---|---|---|---|---|---|---|---|

Accuracy error (%) | Precision error (%) | Accuracy error (%) | Precision error (%) | ||||

CrW | PhD—SNR 24 dB | 0.69 | 6.23 | 11.67 | 14.39 | 2.78 | 0.73 |

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 |

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 |

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 |

TofF—SNR 15 dB | 1.25 | 3.83 | 1.42 | 1.68 | 11.67 | 1.25 |

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. $R2080$ 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\u2009\u2009\mu 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.

**F16 :**

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$.

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 |

Gel 10% | Gel 15% | CNR | $R2080$ (mm) | Time (s) | |||
---|---|---|---|---|---|---|---|

Accuracy error (%) | Precision error (%) | Accuracy error (%) | Precision error (%) | ||||

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 |

**F17 :**

(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: $cshear\u200910%=3.11\xb10.05\u2009\u2009m/s$, and $cshear\u200915%=4.78\xb10.06\u2009\u2009m/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\xb10.82\u2009\u2009m/s$ and $7.47\xb10.52\u2009\u2009m/s$, respectively. The CNR was calculated to be 1.85.

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.^{17}

Lateral resolution is found to be similar ($R2080<1\u2009\u2009mm$) for CrW, SWP, and TBP; as well as ($R2080<1.4\u2009\u2009mm$) 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.,^{36} 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.^{36} 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.^{37} Here, the estimated shear modulus of the left side of a chicken tibialis anterior muscle of $25.3\xb12.2\u2009\u2009kPa$ can be translated to a shear wave speed of $5.03\xb10.22\u2009\u2009m/s$ using Eq. (6), which is similar to our results for normal muscle ($5.71\xb10.82\u2009\u2009m/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.,^{38} and Song et al.^{39} In addition, holographic techniques^{40}^{,}^{41} 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^{42} and human cornea.^{43} 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^{44} offers the possibility to resolve depth-layered media using Rayleigh wave propagation by solving a constrained inverse problem.

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

**Sandrin L.**et al., “Transient elastography: a new noninvasive method for assessment of hepatic fibrosis,” Ultrasound Med. Biol.. 29, (12 ), 1705 –1713 (2003). 0301-5629

**CrossRef**

**Schmitt C.**et al., “Noninvasive vascular elastography: toward a complementary characterization tool of atherosclerosis in carotid arteries,” Ultrasound Med. Biol.. 33, (12 ), 1841 –1858 (2007). 0301-5629

**CrossRef**

**Rogowska J.**et al., “Optical coherence tomographic elastography technique for measuring deformation and strain of atherosclerotic tissues,” Heart. 90, (5 ), 556 –562 (2004).

**CrossRef**

**Donald B. P.**et al., “Visualization and quantification of breast cancer biomechanical properties with magnetic resonance elastography,” Phys. Med. Biol.. 45, (6 ), 1591 –1610 (2000). 0031-9155

**CrossRef**

**Detorakis E. T.,**and

**Pallikaris I. G.**, “Ocular rigidity: biomechanical role, in vivo measurements and clinical significance,” Clin. Exp. Ophthalmol.. 41, (1 ), 73 –81 (2013).

**CrossRef**

**Parker K. J.,**,

**Doyley M. M.,**and

**Rubens D. J.**, “Imaging the elastic properties of tissue: the 20 year perspective,” Phys. Med. Biol.. 56, (1 ), R1 (2011). 0031-9155

**CrossRef**

**Muthupillai R.**et al., “Magnetic resonance elastography by direct visualization of propagating acoustic strain waves,” Science. 269, (5232 ), 1854 –1857 (1995). 0036-8075

**CrossRef**

**Drexler W.**et al., “Optical coherence tomography today: speed, contrast, and multimodality,” J. Biomed. Opt.. 19, (7 ), 071412 (2014). 1083-3668

**CrossRef**

**Mulligan J. A.**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). 1077-260X

**CrossRef**

**Schmitt J. M.**, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express. 3, (6 ), 199 –211 (1998). 1094-4087

**CrossRef**

**Wang R. K.,**,

**Ma Z.,**and

**Kirkpatrick S. J.**, “Tissue Doppler optical coherence elastography for real time strain rate and strain mapping of soft tissue,” Appl. Phys. Lett.. 89, (14 ), 144103 (2006). 0003-6951

**CrossRef**

**Li C.**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). 1742-5689

**CrossRef**

**Manapuram R.**et al., “Estimation of shear wave velocity in gelatin phantoms utilizing PhS-SSOCT,” Laser Phys.. 22, (9 ), 1439 –1444 (2012). 1054-660X

**CrossRef**

**Song S.,**,

**Huang Z.,**and

**Wang R. K.**, “Tracking mechanical wave propagation within tissue using phase-sensitive optical coherence tomography: motion artifact and its compensation,” J. Biomed. Opt.. 18, (12 ), 121505 (2013). 1083-3668

**CrossRef**

**Li C.**et al., “Full skin quantitative optical coherence elastography achieved by combining vibration and surface acoustic wave methods,” Proc. SPIE. 9322, , 93220O (2015). 0277-786X

**CrossRef**

**Nahas A.**et al., “From supersonic shear wave imaging to full-field optical coherence shear wave elastography,” J. Biomed. Opt.. 18, (12 ), 121514 (2013). 1083-3668

**CrossRef**

**Song S.**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). 0146-9592

**CrossRef**

**Wang S.,**and

**Larin K. V.**, “Noncontact depth-resolved micro-scale optical coherence elastography of the cornea,” Biomed. Opt. Express. 5, (11 ), 3807 –3821 (2014). 2156-7085

**CrossRef**

**Li J.**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). 1083-3668

**CrossRef**

**Wang S.,**and

**Larin K. V.**, “Shear wave imaging optical coherence tomography (SWI-OCT) for ocular tissue biomechanics,” Opt. Lett.. 39, (1 ), 41 –44 (2014). 0146-9592

**CrossRef**

**Ahmad A.**et al., “Magnetomotive optical coherence elastography using magnetic particles to induce mechanical waves,” Biomed. Opt. Express. 5, (7 ), 2349 –2361 (2014). 2156-7085

**CrossRef**

**Pinton G. F.,**,

**Dahl J. J.,**and

**Trahey G. E.**, “Rapid tracking of small displacements with ultrasound,” IEEE Trans. Ultrasonics Ferroelectrics Freq. Control. 53, (6 ), 1103 –1117 (2006).

**CrossRef**

**Loupas T.,**,

**Peterson R. B.,**and

**Gill R. W.**, “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).

**CrossRef**

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

**Elegbe E. C.,**and

**McAleavey S. A.**, “Single tracking location methods suppress speckle noise in shear wave velocity estimation,” Ultrasonic Imaging. 35, (2 ), 109 –125 (2013). 0161-7346

**CrossRef**

**Glozman T.,**and

**Azhari H.**, “A method for characterization of tissue elastic properties combining ultrasonic computed tomography with elastography,” J. Ultrasound Med.. 29, (3 ), 387 –398 (2010). 0278-4297

**CrossRef**

**Zvietcovich F.**et al., “Experimental classification of surface waves in optical coherence elastography,” Proc. SPIE. 9710, , 97100Z (2016). 0277-786X

**CrossRef**

**Wu Z.**et al., “Sonoelastographic imaging of interference patterns for estimation of the shear velocity of homogeneous biomaterials,” Phys. Med. Biol.. 49, (6 ), 911 –922 (2004). 0031-9155

**CrossRef**

**Hah Z.**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). 0301-5629

**CrossRef**

**Hoyt K.,**,

**Castaneda B.,**and

**Parker K. J.**, “Two-dimensional sonoelastographic shear velocity imaging,” Ultrasound Med. Biol.. 34, (2 ), 276 –288 (2008). 0301-5629

**CrossRef**

**Schmidt A.,**and

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

**Yao J.**et al., “Angular scan optical coherence tomography imaging and metrology of spherical gradient refractive index preforms,” Opt. Express. 23, (5 ), 6428 –6443 (2015). 1094-4087

**CrossRef**

**Meemon P.,**,

**Lee K.,**and

**Rolland J.**, “Doppler imaging with dual-detection full-range frequency domain optical coherence tomography,” Biomed. Opt. Express. 1, (2 ), 537 –552 (2010). 2156-7085

**CrossRef**

**Meemon P.**et al., “Crawling wave optical coherence elastography,” Opt. Lett.. 41, (5 ), 847 –850 (2016). 0146-9592

**CrossRef**

**Rouze N. C.**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).

**CrossRef**

**Ormachea J.**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).

**CrossRef**

**Koo T. K.**et al., “Relationship between shear elastic modulus and passive muscle force: an ex-vivo study,” J. Biomechan.. 46, (12 ), 2053 –2059 (2013).

**CrossRef**

**Singh M.**et al., “Phase-sensitive optical coherence elastography at 1.5 million A-Lines per second,” Opt. Lett.. 40, (11 ), 2588 –2591 (2015). 0146-9592

**CrossRef**

**Song S.**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). 0003-6951

**CrossRef**

**Liu C. H.**et al., “Non-contact single shot elastography using line field low coherence holography,” Biomed. Opt. Express. 7, (8 ), 3021 –3031 (2016). 2156-7085

**CrossRef**

**Li S.**et al., “Toward soft-tissue elastography using digital holography to monitor surface acoustic waves,” J. Biomed. Opt.. 16, (11 ), 116005 (2011). 1083-3668

**CrossRef**

**Nightingale K. R.**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).

**CrossRef**

**Lombardo G.**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). 1932-6203

**CrossRef**

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", *J. Biomed. Opt*. 22(3), 035010 (Mar 30, 2017). ; http://dx.doi.org/10.1117/1.JBO.22.3.035010

## Figures

**F2 :**

(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)].

**F3 :**

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

**F4 :**

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.

**F6 :**

(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.