Research Papers: General

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

[+] Author Affiliations
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
History: Received November 18, 2016; Accepted March 15, 2017
Text Size: A A A

Open Access Open Access

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 (10 to 100 mm) with resolutions on the order of 10  mm.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,1215 acoustic radiation force,16,17 air-puff excitation,1820 and magnetic transducers21 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 Larin18 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 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,t1)ϕ(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 Display Formula

uz(z)=Δϕ(z)λ04πn,(1)
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 uz(x,z), in which a mechanical wave (i.e., shear wave) propagating in the medium can be tracked.

Graphic Jump Location
Fig. 1
F1 :

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

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

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 Display Formula

I(z,t0)=|I(z,t0)|ejϕ(z,t0),(2)
Display Formula
I(z,t1)=|I(z,t1)|ejϕ(z,t1),(3)
where |I| represents the magnitude, ϕ is the phase of the signal at two consecutive instants t0 and t1, and t0<t1. A direct way to estimate the phase difference Δϕ(z)=  ϕ(z,t1)ϕ(z,t0) can be achieved by Display Formula
Δϕ(z)=arg[I(z,t1)I(z,t0)*],(4)
Display Formula
Δϕ(z)=arctan{Im[I(z,t1)I(z,t0)*]Re[I(z,t1)I(z,t0)*]},(5)
where I* is the complex conjugated version of I, and Im{·} and Re{·} 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.,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, Δϕ(z) can be represented as Display Formula
Δϕ(z)=arctan{m=0M1[Q(zm,t0)II(zm,t1)II(zm,t0)Q(zm,t1)]m=0M1[II(zm,t0)II(zm,t1)+Q(zm,t0)Q(zm,t1)]}1+arctan{m=0M2n=01[Q(zm,tn)II(zm+1,tn)II(zm,tn)Q(zm+1,tn)]m=0M2n=01[II(zm,tn)II(zm+1,tn)+Q(zm,tn)Q(zm+1,tn)]}/2π,(6)
where zm is the element m’th along the z-axis, II(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)].

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.

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 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 Display Formula

W(x,z;t)=Aeα(x+D2)ej[k(x+D2)ωt],(7)
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.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 Display Formula
cshear=ωk=λf,(8)
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 as6Display Formula
cshear=Gρ.(9)

Graphic Jump Location
Fig. 2
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)].

In general, biological tissues are nearly incompressible and the Poisson’s ratio can be approximated to ν=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 λ or, equivalently, by estimating the wave number k=2π/λ. In addition, the shear wave speed can also be estimated by using the kinematic relationship Display Formula

cshear=ΔdshearΔtshear,(10)
where Δdshear is the propagated distance of a transient shear wave tone during time Δtshear. This estimation technique is also called the time-of-flight method and is widely used in ultrasound (US) and OCT elastography.25 Shear wave velocity in tissue typically ranges from 0.1 to 10  m/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 ν=0.5 is approximately given by27Display Formula

cshear1.05*cRayleigh.(11)
For further details about Rayleigh waves in OCE, we refer the reader to our previous study by Zvietcovich et al.27

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

Wleft(x,t)=Aeα(x+D2)ej[k(x+D2)ωt+φl],(12)
Display Formula
Wright(x,t)=Aeα(D2x)ej[(k+Δk)(D2x)(ω+Δω)t+φr],(13)
where k and k+Δk are the wave numbers of Wleft and Wright, 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=±D/2, respectively. According to Wu et al.,28 assuming no reflections in the medium, the square amplitude of the interference U(x,t)=Wleft+Wright is given by Display Formula
|U(x,t)|2=2A2eαD{cosh(2αx)+cos[Δωt+2(k+Δk2)xφ0]}.(14)
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) pattern28 [Fig. 2(c)]. The phase term φ0=Δk2D+Δφ is constant and Δφ=φrφl. The CrW speed cCrW is related to the shear wave speed cshear by Display Formula
cCrW=cshearΔω2ω+Δω,(15)
where cshear=ω/k, and cCrW=Δω/(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 cshear 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, cCrW=0, and the CrW pattern becomes a standing wave (SW). In this situation, Display Formula
|U(x,t)|2=2A2eαD[cosh(2αx)+cos(2  kxφ0)],(16)
and the shear wave speed estimation reduces to the estimation of k in the SW pattern.

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

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.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 Δω in order to remove the nonharmonic cosh term. Then, assuming Δωω, it is also true that Δkk, and Display Formula

|U(x,t)|2=B{cos[Δωt+ϕ(x)]},(17)
where B=2A2eαD is the amplitude and ϕ(x)=2kxφ0 is the phase signal for each lateral position. The estimated version of the phase ϕ˜(x) can be computed by means of harmonic fitting or by using the FT approach. The estimated shear wave speed c˜shear(x) can be calculated using Display Formula
c˜shear(x)=ωk˜=4πfϕ˜(x)  ,(18)
where ϕ˜(x) is the spatial derivative of ϕ˜(x), and ω=2πf. The derivative operator is vulnerable to noise corruption; then a smoothing procedure can be applied to ϕ˜(x) at the expense of losing lateral resolution. c˜shear(x) is computed for each depth z in order to create a 2-D shear wave velocity map.

Hoyt’s estimator

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 Display Formula

R(n,1)=1L1m=0L2sA*(n+m)sA(n+m+1).(19)

Therefore, the spatial wave number can be estimated as Display Formula

k(xn)=1Δxarctan{Im[R(n,1)]Re[R(n,1)]},(20)
where xn=nΔx, for n=0,1,2,,N1L and Δx is the lateral sampling resolution. Using Eqs. (20) and (18), a 2-D shear-wave-speed map (SWSM) can be computed.

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

cpeak(xn)=[i=0w1pi(n)qi(n)]wp(n)¯  q(n)¯[i=0w1qi(n)2]wq(n)¯2,(21)
where q(n)¯=1wi=0w1qi(n), p(n)¯=1wi=0w1pi(n), xn=nΔx, tn=nΔt for n=0,1,2,,Nw, 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.

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.

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 1.35, which is close to that of real tissue. The shape of the phantom is a rectangular cuboid of dimensions 15×7×3  cm3 [Fig. 3(a)].

The Young’s modulus of each side of the phantom was measured by taking three cylindrical samples (n=3) 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.31 Finally, using Eq. (9) for a material density of 1  g/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.

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/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  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.32 The imaging depth of the system was measured to be 5 mm in air (10  dB sensitive fall-off). The optical lateral resolution was 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.33 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 (uz,min) by plugging Δϕerr in Eq. (1). We found a displacement of uz,min=0.358  nm. Finally, the maximum axial displacement supported by the system without unwrapping the phase-shift signal (Δϕmax=π) is uz,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 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 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 uz,max=0.24  μm without unwrapping the Doppler phase signal. Since waves travel 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.

Graphic Jump Location
Fig. 3
F3 :

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

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

Graphic Jump Location
Fig. 4
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.

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.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  Hz and fright=400.4  Hz, respectively, yielding a Δf=0.4  Hz or Δω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’ approach23 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.

Graphic Jump Location
Fig. 5
F5 :

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

Graphic Jump Location
Fig. 6
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.

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 x0, M depth A-lines are acquired at a rate of 20 kHz (time resolution ΔM=50  μs) as shown in Fig. 7.

Graphic Jump Location
Fig. 7
F7 :

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 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Δt, x=nΔx, for n=0,1,2,,N1, 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 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, Δt=20  ms. Equation (14) can be rewritten as a function of Δx and Δt as Display Formula

|U(n)|2=2A2eαD{cosh(2αnΔx)+cos[ΔωnΔt+2(k+Δk2)nΔxφ0]}.(22)

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 Display Formula

|U(n)|2=B(cos{[ΔωΔt+2(k+Δk2)Δx]nφ0}),(23)
where B=2A2eαD. Then the spatial discrete wave number of the signal is ΔxkSCrW=ΔωΔt+2(k+Δk2)Δx, where kSCrW is the spatial continuous wave number, and Display Formula
kSCrW=ΔωΔtΔx+2(k+Δk2).(24)
Equation (24) enables the recovery of k (needed for estimating the local shear wave speed) based on kSCrW. Moreover, if wave number varies along the lateral axis [k(x)], then kSCrW(x) can be found using two estimators: local fitting and Hoyt’s algorithm.30Figure 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.

Graphic Jump Location
Fig. 8
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.

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 Display Formula

|U(n)|2=2A2eαD[cosh(2αnΔx)+cos(2knΔxφ0)].(25)
Using the appropriate filter, the DC bias can be suppressed, so that Display Formula
|U(n)|2=B[cos([2kΔx]nφ0)].(26)
Then the SW number kSW=2k can be recovered by applying Hoyt’s method30 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).

Shear wave propagation method

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  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 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ΔM, and x=nΔx, for m=0,1,2,,M1, and n=0,1,2,,N1. Δx and Δ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 Δg=5  ms needs to be added to t. Then t=mΔM+nΔg, and Eq. (7) becomes Display Formula

W(m,n)=Aeα(nΔx+D2)ej[k(nΔx+D2)ωmΔM]ej(ωnΔg),(27)
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: Display Formula
β(n)=β(n+1).(28)
Then, ωΔg=2πl is needed and lN+. 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×103)=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 ϕ(xn)=knΔx can be retrieved from Eq. (27) [Fig. 9(b)], and the 2-D shear wave speed may be computed.

Graphic Jump Location
Fig. 9
F9 :

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

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.

Graphic Jump Location
Fig. 10
F10 :

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)×1005(depth axis) elements. We rearranged the scanning into a three-dimensional (3-D) volume of 1000×1005×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).

Graphic Jump Location
Fig. 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).

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 cs1=3.11  m/s, and cs2=4.78  m/s, respectively. The size of the 2-D grid is 78  mm(600elements)×43  mm(300elements) 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.

Graphic Jump Location
Fig. 12
F12 :

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.

Graphic Jump Location
Fig. 13
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 μm related to the particle displacement.

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.35 Then the speed profile is fitted to the function Display Formula

c(x)=c10%+(c15%c10%)11+ex0xτ,(29)
where c10% and c15% are the wave velocities in the 10% and 15% regions, respectively, x0 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 R2080, defined as the distance from 20% to the 80% of the transition.35 Then the lateral resolution is equivalent to Display Formula
R2080=2τln(4).(30)

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 Display Formula

CNR=|μ15%μ10%|σ15%2+σ10%2,(31)
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×4  mm2 for the experimental results and 8×8  mm2 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 Display Formula

Ea=|μ10/15%μgt|μgt100%.(32)
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 Display Formula
Ep=σ10/15%μ10/15%100%.(33)

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 Display Formula

σ+ηE1Dtασ=E0ϵ+ηE0+E1E1Dtαϵ,(34)
where σ is the stress, ϵ is the strain, Dtα is the fractional derivative operator, and the four parameters of the model: E0 (left spring), E1 (right spring), and η, α (spring-pot) as described in Schmidt et al.31 Once the four parameters are found using the optimization approach, Eq. (34) can be transformed in the frequency domain to Display Formula
σ˜=E0+(i2πf)αηE0+E1E11+(i2πf)αηE1ϵ˜=E*ϵ˜,(35)
where f is the frequency, σ˜ and ϵ˜ 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*=3G* in Eq. (9).

Graphic Jump Location
Fig. 14
F14 :

(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±0.05  m/s and 4.78±0.06  m/s, respectively.

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

Graphic Jump Location
Fig. 15
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.

Table Grahic Jump Location
Table 1Accuracy 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: cshear10%=3.11±0.05  m/s and cshear15%=4.78±0.06  m/s.
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. 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  μ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.

Graphic Jump Location
Fig. 16
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.

Table Grahic Jump Location
Table 2Mean and standard deviation of shear wave velocity estimation for all methods in both phantom gelatin concentrations. Ground truth values: cshear10%=3.11±0.05  m/s and cshear15%=4.78±0.06  m/s.
Table Grahic Jump Location
Table 3Accuracy error and precision of shear wave speed estimation, CNR, resolution (R2080), and acquisition time for all methods extracted from experimental results. Ground truth values: cshear10%=3.11±0.05  m/s and cshear15%=4.78±0.06  m/s.
Graphic Jump Location
Fig. 17
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: cshear10%=3.11±0.05  m/s, and cshear15%=4.78±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±0.82  m/s and 7.47±0.52  m/s, respectively. The CNR was calculated to be 1.85.

Graphic Jump Location
Fig. 18
F18 :

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.

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  mm) for CrW, SWP, and TBP; as well as (R2080<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.,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±2.2  kPa can be translated to a shear wave speed of 5.03±0.22  m/s using Eq. (6), which is similar to our results for normal muscle (5.71±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.,38 and Song et al.39 In addition, holographic techniques40,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 liver42 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 Oldenburg44 offers the possibility to resolve depth-layered media using Rayleigh wave propagation by solving a constrained inverse problem.

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

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

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
Mohan  K., and Oldenburg  A., “Elastography of soft materials and tissues by holographic imaging of surface acoustic waves,” Opt. Express. 20, , 18887 –18897 (2012). 1094-4087 CrossRef
© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

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


Figures

Graphic Jump Location
Fig. 2
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)].

Graphic Jump Location
Fig. 1
F1 :

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

Graphic Jump Location
Fig. 7
F7 :

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.

Graphic Jump Location
Fig. 8
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.

Graphic Jump Location
Fig. 10
F10 :

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