Accuracy of Common Motion Estimators in Wave-Based Optical Coherence Elastography

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.


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

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

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: Here, -%&'( 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 . %)*+,-& 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.

Motion Estimation Algorithms
The Kasai autocorrelation algorithm can be used to calculate the average phase change over a given axial range, Δϕ 1111 : 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 i th A-scan at a given depth [10,11].  The Loupas autocorrelator can similarly be used to find Δϕ [1,10]: where 2 . is the m th 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 where 9 ! is the first frame, 9 / is the second frame, and 5 = 1 … " 1 and 7 = 1 … " 2 are the horizontal and vertical indices of the processing window, respectively. Then, 3 is averaged laterally across a window of the desired length, giving the vector B(j): where Φ is the phase variation. After this, complexvalued quantities c(j) are formed via: and averaged across the vertical index as: where C is the complex vector. The angle ΔΦ 1111 is the desired quantity and is the averaged inter-pixel phasevariation across the processing window.

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.  the 50 depth-wise pixels. The root-mean-square (RMS) error of the individual profiles as compared to their noiseadded 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. 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.

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

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   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, . 78 , is related to the OCT SNR by 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: 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., < ±π), 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 Δ2π , 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 Δ2π 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 crosscorrelation 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/quasistatic 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.

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 wavebased elastography. Our results show that there is a lower J of Biomedical Photonics & Eng 7(4) DD Mm 2021 © J-BPE 040303-7 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.