CC BY
13
Accuracy of Common Motion Estimators in Wave-Based Optical Coherence Elastography
Justin R. Rippy1*, Manmohan Singh1, Salavat Aglyamov2, and Kirill V. Larin1
1 Department of Biomedical Engineering, University of Houston, 3517 Cullen Blvd., Houston, TX 77204, USA
2 Department of Mechanical Engineering, University of Houston, Houston, TX 77204, USA
* e-mail: [email protected]
Abstract. Motion estimators are commonly used in shear wave optical coherence elastography to compute small displacements. This work focuses on comparing three motion estimators: Kasai, Loupas, and the vector method. Our results show that the vector method is superior in low-SNR, low-amplitude situations for group velocity calculation and profile reconstruction. © 2021 Journal of Biomedical Photonics & Engineering.
Keywords: motion estimator; shear wave elastography; signal to noise ratio; Loupas; Kasai; vector method.
Paper #3450 received 4 Aug 2021; revised manuscript received 7 Sep 2021; accepted for publication 9 Sep 2021; published online 21 Sep 2021. doi: 10.18287/JBPE21.07.040303.
1 Introduction
In traditional shear wave elastography, mechanical waves are generated in a sample and measured by various imaging modalities [1]. One modality that is often used for this purpose is optical coherence tomography, which utilizes backscattered light to obtain images with micrometer resolution [2]. By determining the speed of the mechanical wave traveling through a material, mechanical properties can be estimated, which can be a powerful tool for early disease detection, disease progression monitoring, or therapy guidance [3]. Small displacements on the order of nanometers to micrometers are typically induced in optical coherence elastography (OCE), and the ability to evaluate these displacements determines the quality of speed estimation regardless of the speed calculation method. Motion estimators are often used for this task, which computes the motion between temporal and/or spatial samples utilizing the raw data acquired from the imaging system. Three of the most often used motion estimators are the Kasai [4] and Loupas [5] autocorrelation estimators and the recently introduced vector method [6]. The Kasai autocorrelator, which is also known as the 1D autocorrelator, was used extensively in initial elastography publications due to its computational simplicity [5]. Later, Loupas was able to derive a superior estimator known as the 2D autocorrelator, which deals with some of the limitations of the Kasai method (which are discussed in Ref. [5]). The most recent of the three, the vector method tracks subpixel displacements between scans to ensure that the same scatterers are compared when determining displacement [6, 9]. This method has been shown to be
more accurate and more noise-tolerant than other methods but only for static/quasi-static elastography [6, 9]. However, its application to wave-based elastography has yet to be demonstrated. Though these three methods have been and continue to be used throughout elastography, there are no studies showing which motion estimator is best for a given set of circumstances.
In this work, we performed simulations and collected experimental data to determine what effect varying the amplitude and signal-to-noise ratio (SNR) had upon the measured group velocity. We used normalized fitted profiles to create ideal displacement profiles with infinite SNR to be used as a baseline for comparison. White Gaussian noise was added at various levels to match desired SNR values. Profiles were reconstructed using the three methods mentioned earlier and using only the raw data with and without averaging. Following this, the group velocity of the mechanical wave was calculated, and the results of each method were compared to the normalized fitted profiles. We show that there is an apparent lower bound on the SNR required to accurately obtain the wave speed with a given speed calculation method. We also show that the vector method is more resistant to the effects of noise than the Loupas or Kasai motion estimators but at the cost of a slight loss of accuracy.
2 Materials and Methods
2.1 Data acquisition
Elastography data was acquired from several gelatin phantom concentrations (8%, 12%, and 16%) with an optical coherence elastography imaging system. The system and acquisition procedure is described in a previous publication [7]. The OCE system contains a broadband swept-source laser (HSL 2000, Santec Corp., Hackensack, NJ, USA) with a central wavelength of 1310 nm, the axial resolution of 11 ^m in air, transverse resolution of 16 ^m, and a scan rate of 30 kHz. A focused micro air pulse [8] was synchronized with the frame trigger [9] to give an effective temporal resolution of 30 kHz. This air pulse generated a wave on the surface of the phantom, which was then recorded at 251 points in a line over ~6.3 mm. The excitation was performed in the center of this line, causing wave motion to the left and right. The OCE data was obtained as the raw spectral interference pattern, which was then converted to linear-in-wavenumber spectral data by linear interpolation based on calibration. The raw complex OCT data was then obtained by fast Fourier transform (FFT) of the linear-in-wavenumber data. In the following simulations, profiles were generated from this experimental OCT data at each concentration mentioned above.
2.2 Simulations
From the raw complex OCT data, the top surface of each phantom was determined using a surface tracking algorithm, and displacement profiles were taken for 50 pixels (~200 ^m) below the surface. Stationary phase profiles were obtained from the experimental OCT phase data (i.e., angle of the raw complex OCT data). A Gaussian fit was performed on the stationary phase profiles to generate ideal noiseless profiles used in the following simulations. Following this, only the ideal noiseless profiles generated from the Gaussian fit were used. The width of the fit was fixed to eliminate any effect of dispersion on the results. For these simulations, the amplitudes and SNRs of the displacement profiles were modified before motion estimation. The displacement profiles were first normalized and multiplied by the desired amplitude to give profiles with infinite SNR but definite amplitude. An example profile corresponding to an amplitude of -0.1 is shown in Fig. 1.
To modify the SNR, white Gaussian noise was added to each profile using until the desired SNR was obtained, where the SNR was given by:
SNR = 20 log10
\$profile)
(1)
Here, upeak is the average of the 3 central samples at the known maximum peak location obtained from the fitted profiles (e.g., the minimum in Fig. 1), and .profue is the standard deviation of the profile after subtraction of the fitted profile.
Several example profiles with varying SNRs are shown in Fig. 2. Twelve individual amplitudes and SNR values were combined to study their effects on group velocity calculation. The amplitudes and SNRs investigated are listed in Table 1. From these noise-added profiles, the various motion estimators were applied: no motion estimator was applied (i.e., only the angle of the complex data was utilized), smoothing the noise-added profiles in-depth and time, Kasai estimator, Loupas estimator, and vector method. For all techniques, a kernel size of 5 pixels was utilized. Since some techniques, such as the Kasai and Loupas estimators, provide the particle velocity, i.e., the time derivative of the displacement, the cumulative sum was applied to obtain the displacement from the particle velocity.
Normalized fitted phase profile
8 10 12 Time (ms)
Fig. 1 A normalized fitted phase profile with amplitude 1 is shown here.
Table 1 Amplitude and SNR values used in the simulations.
Amplitude (radians)
0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 5, 7.5, 10, 25, 50, 100
SNR (dB)
0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 5, 7.5, 10, 25, 50, 100
3 Motion Estimation Algorithms
The Kasai autocorrelation algorithm can be used to calculate the average phase change over a given axial range, Acjx
Ac# :
-tan
_i iï,?![Q(iV(i-l)-l(i)Q(i-l)]\ (2) Z#="i1[/(i)/(i-i) + Q(i)Q(i-i)]
where Q and I are the imaginary and real parts of the complex OCT signal, respectively; N is the number of averaged A-scans, and i is the signal of the ith A-scan at a given depth [10, 11].
SNR: 1.5 dB
SNR: 2.5 dB
200
300 400 Time (ms)
SNR:10 dB
500
200
300 400 Time (ms)
500
200 300 400 Time (ms)
SNR:100 dB
500
200
300 400 Time (ms)
500
Fig. 2 Example profiles at various SNR values. Each profile started with an amplitude of 1.
The Loupas autocorrelator can similarly be used to find A< [1, 10]:
A^(z) =
tan
"Z!&1o[Q(zm,to)I(zm,t1)-I(zm,t0)Q(zm,t1)]\ Z10[I(zm,t0)I(zm,t1)+Q(zm,t0)Q(zm,t1)]J
-1 <Tm=l\
vzm=ii
Zm=iZ'.= o[>(+m.tn)l(Zm+1.tn) + Q(zm.tn)Q(Zm+1.tn)]J
1 +-
(3)
26
where z. is the mth element along the z-axis, M is the size of the window in the axial direction, N is the size of the window in time, and Q and I are the imaginary and real parts of the complex OCT signal, respectively.
In the vector method, there are four steps [6]. First, the complex data from two temporally consecutive OCT frames are multiplied together according to
b(m,j) = a1(m,j) • a/(m,j),
(4)
where a1 is the first frame, a2 is the second frame, and m=l...Nx and j = l...Nz are the horizontal and vertical indices of the processing window, respectively. Then, b is averaged laterally across a window of the desired length, giving the vector B(j):
b0) = B(j) exp[i ■ 00')] = b(m,j), (5)
where 0 is the phase variation. After this, complex-valued quantities c(j) are formed via:
c(j) = b{] + I) ■ b{])*, (6)
and averaged across the vertical index as:
c = r^!1 c(j) = Cexp[i
AO],
(7)
where C is the complex vector. The angle AO is the desired quantity and is the averaged inter-pixel phasevariation across the processing window.
4 Group Velocity Calculations
A cross-correlation-based algorithm was utilized to quantify the elastic wave group velocity [12]. The temporal phase profiles along the propagation path for a given depth were cross-correlated with a reference signal near the excitation to obtain the wave propagation lags. The least-squares regression fitting of the wave propagation lags with the propagation distances resulted in the wave speed. The speed was calculated for the various motion estimators: noiseless fitted profiles, noise-added data, averaged noise-added profiles, Kasai estimated motion, Loupas reconstructed profiles and vector method data.
5 Comparison with Noiseless Profiles
Following the wave propagation simulations, the results of each method were compared with the results of the noiseless profiles. For the 8%, 12%, and 16% phantoms, the group velocities of the noiseless fitted profiles were 2.518 m/s, 3.139 m/s, and 8.085 m/s, respectively. The mean and standard deviation of the group velocities for each method were calculated and used to calculate percent error for each method, with the group velocity generated from the noiseless profiles serving as the true value. The mean and standard deviation were taken for
Fig. 3 Comparison of the RMS error of the three motion estimators from 8%, 12%, and 16% gelatin phantom data.
the 50 depth-wise pixels. The root-mean-square (RMS) error of the individual profiles as compared to their noise-added counterparts, which was the given motion in the samples, was also calculated to estimate the accuracy of each method to reconstruct the given motion. Finally, a summary of the SNR improvement was also quantified, where the SNR is given by Eq. 4.
6 Results
Fig. 3 shows the RMS error between the reconstructed phase profiles and the noiseless profiles for each method. In these figures, RMS error is capped at 20% to provide better contrast. The vector method has a larger area of RMS error less than 5% (indicated by the blue region) than either of the other two methods. Kasai outperforms Loupas here as well. In the low-SNR region, the vector method largely performs worse than either method except at very low SNR and amplitude.
Group velocity results for the three motion estimation methods are shown in Fig 4. The percent error is capped at 100% in these figures to provide better contrast for the lower values. The average percent error across all three phantom concentrations at SNRs greater than 50 was 1.94 ± 1.56%, 2.07 ± 1.44%, and 2.30 ± 2.10% for the Kasai, Loupas, and vector methods, respectively. For an
SNR of 25, only the amplitudes less than 10 rad had percent error values less than 10%, at 3.09 ± 2.29% for Kasai, 2.3 ± 1.62% for Loupas, and 2.39 ± 2.22% for the vector method. However, for SNR values lower than 10 dB all methods had difficulty accurately determining wave speed except at extremely low amplitudes (< 1 - 2 rad), with percent errors sometimes exceeding 900%. The vector method accurately determined the group velocity in the areas where low SNR meets with low amplitudes in all cases, unlike either Kasai or Loupas. This can be seen in the lower-left corner of Fig. 4 in plots (c), (f), and (i).
7 Discussion
All three methods allowed for highly accurate (< 5% error) group velocity calculations regardless of amplitude above an OCE SNR of ~50 dB. However, at SNRs above 50 dB, motion estimation algorithms are not necessary, as the group velocities can be calculated directly from the raw phase data, as shown in Fig 5. This may not be the case for samples with more complicated excitation patterns (e.g., chirps) or in cases with irregular noise. More complicated excitation patterns, especially those with multiple excitation waveforms, may present difficulty in retrieving group velocity due to errors in
J
20 40 60 SNR(dB)
Kasai, 12%
Percent Error 100.0
- 90.00
- 80 00
- 70.00 60 00
- 50.00 40 00 30 00 20.00 10 00 0.000
Percent Error
Loupas, 8%
J
Percent Error C ) H- 100.0 1001
Vector method, 8%
20 40 60 60 SNR (dB)
Kasai, 16%
■ - 100.0 90.00 80.00 70.00 TD 100 80-
- 60.00 S. 60
50.00 TJ £
- 40.00 Q. E < 40-
30.00
_ - 20.00 20-
10.00
M - 0.000
100
Loupas, 12%
y
^ 90.00 - 80 00
- 70.00
- 60.00
- 50.00 40.00
- 30.00 20.00
^ 10.00 I 0.000
Percent Error
n 100.0 90.00 80.00 70.00 - 60.00
Percent Error
H 100.0 '100
H- 90.00 ™ 80.00 80
- 70 .00 Ç
60.00 60 «
I - 50.00 "§ 40.00
- 30.00
J 20.00 I- 10.00 I- 0.000
Loupas, 16%
20 40 60 80 SNR (dB)
J
-40.00 B. 40-
Percent Error ^m 100.0 90.00 80.00
■ 70.00 60.00
■ 50.00 40.00
■ 30.00 20.00 10.00 0.000
Vector method, 12%
Vector method, 16%
Percent Error ioo o 90.00 - 80.00 - 70.00 60.00 50.00 40.00 30.00 20.00 10.00 0.000
Percent Error
■ 100.0
90.00 - 80.00 70.00 - 60.00 I 50.00 - 40.00 30.00
Î 20.00 10.00 0.000
Percent Error ^m 100.0 ^T 90-00 60.00
- 70.00 60.00 50.00
- 40.00
- 30.00
- 20.00 - 10.00 - 0.000
20 40 60 60 SNR (dB)
Fig. 4 Comparison of the percent error of the three motion estimators from 8%, 12%, and 16% gelatin phantom data. Gray areas indicate that values could not be calculated in that region.
i motion estimator,
J
J
Percent Error 100.0 ■ 90.00 80.00 70 00 60.00 50.00 - 40.00 30.00 20 00 10.00 0.000
Fig. 5 Comparison of the percent error using (a) no motion estimator and (b) after averaging in-depth and time.
cross-correlation. Hence, other techniques for velocity calculation can be applied, such as frequency-domain methods. This result indicates that the SNR can be purposefully lowered by adjusting other parameters such as detector speed without any penalty in velocity error, provided that the SNR remains greater than 50 dB.
In this paper, we discuss the OCE SNR (i.e., the SNR of the displacement of the elastic wave), not the OCT
SNR. However, Park showed in Ref. [13] that the minimum detectable phase difference, is related to the OCT SNR by
cta$ = (SNR)'-.
(8)
This means that having a high OCT SNR will contribute to having a high OCE SNR.
As shown in Fig. 4, averaging the phase information in-depth and time did not improve group velocity calculation. Though the vector method also involves similar averaging, it averages the complex data. The averaging of the real and imaginary parts of the vectors instead of directly averaging their angles makes the vector method robust and noise-tolerant [6].
There is also a lower bound on SNR required for accurate group velocity retrieval regardless of the method or amplitude at ~20 dB. Above this SNR, at least some amplitudes allow for group velocity retrieval, but below 20 dB, errors increase for all methods. Typical SNRs of modern imaging systems are well above this, but this could be a factor when measuring extremely weak waves, designing low-cost imaging systems with relatively weak lasers in optical coherence tomography and other optical techniques, in portable single-fiber-based probes, etc. The minimum phase difference that can be detected depends on the image SNR and is given by:
= vdfer (8)
This equation assumes that the SNR is >> 1 [13]. At an SNR of ~20 dB and an amplitude of ~20 rad, both the vector method and Loupas allow for reliable group velocity calculations, but Kasai does not. This is likely due to the limitations mentioned previously [5, 10, 14] by virtue of the Kasai autocorrelator being a 1D technique. In the square region from 20 dB to 50 dB and 0 rad to 25 rad, the vector method is more reliable across all wave speeds (i.e., elasticities). This is due to phase unwrapping and the errors that may result from using phase unwrapping algorithms. When the SNR and amplitudes are low, the amplitudes of the profiles remain wrapped (i.e., < ±n), and there are no unwrapping errors, leading to very low error in the final group velocity estimation. As the noise increases and/or the phase jumps increase above A2n , the unwrapping algorithm accumulates errors causing a significant final error in the final group velocity estimation due to large erroneously unwrapped jumps. The focus of this manuscript was to isolate the motion estimators, so various phase unwrapping algorithms were not examined. However, phase unwrapping is still not well solved [15]. It should be given consideration when attempting to reconstruct profiles that are either particularly noisy or that have signals significantly above the A2n threshold. When considering the entirety of the contour map, the vector method typically produces more reasonable group velocities than either Loupas or Kasai techniques, as evident by comparing the regions in Fig. 3. However, this effect is limited to a very small range of SNRs and amplitudes, indicating that either the vector method or the Loupas autocorrelator can be used in the majority of circumstances.
As the wave speed increases, the ability of the estimators to accurately determine velocity also increases. This is likely due to differences in the amount
of peak shifting, also known as jitter, occurring during cross-correlation. Several factors can contribute to these shifts, including signal decorrelation, noise, and finite window length. This causes a jitter in the cross-correlation values computed as the wave propagation lag, which is measured as the standard deviation of the displacement error [10]. As the group velocity slope increases, the error decreases.
The RMS error results show that the vector method reconstructed the motion most accurately, while Kasai and Loupas showed similar performance. This comparison eliminates the effects of cross-correlation and fitting and shows that the vector method is superior to the other methods for profile reconstruction over a higher range of SNRs and amplitudes. Averaging the raw phase profiles did improve the signal, but the vector method was superior because it averages the complex data instead of only the phase data. This accounts for the speckle properties inherent to coherence imaging systems, such as ultrasound imaging and OCT.
The Kasai autocorrelator works by measuring the average phase shift with respect to a central frequency but assumes a constant mean frequency. Because of this assumption, Kasai's algorithm can be inaccurate if there are fluctuations [10]. Additionally, Kasai's algorithm was useful at a time when computing power was limited. Computing power has increased exponentially since the introduction of the algorithm, meaning that more accurate but computationally expensive methods can be used without a significant time penalty.
The Loupas autocorrelator uses information from two dimensions instead of only using axial information for filtering [5]. Hence, the Loupas autocorrelator outperforms the Kasai algorithm. However, in some situations, the denominator of Loupas' algorithm can become small enough that the algorithm yields nonphysical results, and therefore a filter must be used to limit the displacement range to - = to / [10]. The Loupas algorithm is also relatively computationally inexpensive, meaning that results can be obtained quickly. In fact, it was implemented on the GPU for real-time displacement estimation ten years ago [16].
The vector method, which was first published in 2018, maintains low error by using complex data for calculations until the last step and avoiding phase unwrapping [6]. Though it was designed for static/quasi-static elastography, these results show that it can also be utilized for dynamic elastography [6]. Though it is the most computationally intense of the three, numerous computational approaches optimized for vector arithmetic can perform these calculations with minimal time penalties.
8 Conclusions
The primary goal of this work was to look at the relationship between amplitude and SNR when considering motion estimation algorithms for various mechanical wave propagation speeds in impulsive wave-based elastography. Our results show that there is a lower
bound on SNR at approximately 20 dB that is required for accurate group velocity measurement for any of the methods. After an OCE displacement SNR of 50 dB, the percent error remains below 5% regardless of amplitude and motion estimator, which may be useful when trying to balance other parameters that may lower SNR. Additionally, the vector method has greater coverage of lower error values than either Kasai or Loupas methods,
though Loupas is nearly as good. The role of amplitude in profile reconstruction is very small, with SNR being the key parameter, except at very low amplitudes.
Disclosures
All authors declare that there is no conflict of interests in this paper.
References
1. F. Zvietcovich, J. P. Rolland, J. Yao, P. Meemon, and K. J. Parker, "Comparative study of shear wave-based elastography techniques in optical coherence tomography," Journal of Biomedical Optics 22(3), 035010 (2017).
2. C. E. Psomadakis, N. Marghoob, B. Bleicher, and O. Markowitz, "Optical coherence tomography," Clinics in Dermatology, (2021).
3. B. J. Fahey, R. C. Nelson, S. J. Hsu, D. P. Bradway, D. M. Dumont, and G. E. Trahey, "In vivo guidance and assessment of liver radio-frequency ablation with acoustic radiation force elastography," Ultrasound in Medicine & Biology 34(10), 1590-603 (2008).
4. C. Kasai, K. Namekawa, A. Koyano, and R. Omoto, "Real-time two-dimensional blood flow imaging using an autocorrelation technique," IEEE Transactions on Sonics and Ultrasonics 32(3), 458-464 (1985).
5. T. Loupas, R. Peterson, and R. W. Gill, "Experimental evaluation of velocity and power estimation for ultrasound blood flow imaging, by means of a two-dimensional autocorrelation approach," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 42(4), 689-699 (1995).
6. A. Matveyev, L. Matveev, A. Sovetsky, G. Gelikonov, A. Moiseev, and V. Y. Zaitsev, "Vector method for strain estimation in phase-sensitive optical coherence elastography," Laser Physics Letters 15(6), 065603 (2018).
7. J. R. Rippy, M. Singh, S. R. Aglyamov, and K. V. Larin, "Ultrasound Shear Wave Elastography and Transient Optical Coherence Elastography: Side-by-Side Comparison of Repeatability and Accuracy," IEEE Open Journal of Engineering in Medicine and Biology 2, 179-186 (2021).
8. S. Wang, K. V. Larin, J. Li, S. Vantipalli, R. K. Manapuram, S. Aglyamov, S. Emelianov, and M. D. Twa, "A focused air-pulse system for optical-coherence-tomography-based measurements of tissue elasticity," Laser Physics Letters 10(7), 075605 (2013).
9. S. Wang, K. V. Larin, "Shear wave imaging optical coherence tomography (SWI-OCT) for ocular tissue biomechanics," Optics Letters 39(1), 41-44 (2014).
10. G. F. Pinton, J. J. Dahl, and G. E. Trahey, "Rapid tracking of small displacements with ultrasound," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 53(6), 1103-1117 (2006).
11. V. X. Yang, M. L. Gordon, A. Mok, Y. Zhao, Z. Chen, R. S. C. Cobbold, B. C. Wilson, and I. A. Vitkin, "Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation," Optics Communications 208(4-6), 209-214 (2002).
12. S. Wang, A. L. Lopez, Y. Morikawa, G. Tao, J. Li, I. V. Larina, J. F. Martin, and K. V. Larin, "Noncontact quantitative biomechanical characterization of cardiac muscle using shear wave imaging optical coherence tomography," Biomed. Opt. Express 5(7), 1980-1992 (2014).
13. B. H. Park, M. C. Pierce, B. Cense, S.-H. Yun, M. Mujat, G. J. Tearney, B. E. Bouma, and J. F. de Boer, "Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 ^m," Optics Express 13(11), 3931-3944 (2005).
14. D. Morofke, M. C. Kolios, I. A. Vitkin, and V. X. D. Yang, "Wide dynamic range detection of bidirectional flow in Doppler optical coherence tomography using a two-dimensional Kasai estimator," Optics Letters 32(3), 253-255 (2007).
15. K. Wang, Y. Li, Q. Kemao, J. Di, and J. Zhao, "One-step robust deep learning phase unwrapping," Optics Express 27(10), 15100-15115 (2019).
16. S. Rosenzweig, M. Palmeri, and K. Nightingale, "GPU-based real-time small displacement estimation with ultrasound," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 58(2), 399-405 (2011).
