Angular Integral Autocorrelation for Speed Estimation in Shear-Wave Elastography

Abstract

The utilization of a reverberant shear-wave field in shear-wave elastography has emerged as a promising technique for achieving robust shear-wave speed (SWS) estimation. However, many types of estimators cannot accurately measure SWS within such a complicated 3D wave field. This study introduces an advanced autocorrelation estimator based on angular integration known as the angular integral autocorrelation (AIA) approach to address this issue. The AIA approach incorporates all the autocorrelation data from various angles during measurements, resulting in enhanced robustness to both noise and imperfect distributions in SWS estimation. The effectiveness of the AIA estimator for SWS estimation is first validated using a k-Wave simulation of a stiff branching tube in a uniform background. Furthermore, the AIA estimator is applied to ultrasound elastography experiments, magnetic resonance imaging (MRI) experiments, and optical coherence tomography (OCT) studies across a range of different excitation frequencies on tissues and phantoms, including in vivo scans. The results verify the capacity of the AIA approach to enhance the accuracy of SWS estimation and the signal-to-noise ratio (SNR), even within an imperfect reverberant shear-wave field. Compared to simple autocorrelation approaches, the AIA approach can also successfully visualize and define lesions while significantly improving the estimated SWS and SNR in homogeneous background materials and providing improved elastic contrast between structures within the scans. These findings demonstrate the robustness and effectiveness of the AIA approach across a wide range of applications, including ultrasound elastography, magnetic resonance elastography (MRE), and optical coherence elastography (OCE), for accurately identifying the elastic properties of biological tissues in diverse excitation scenarios.

1. Introduction

Shear-wave elastography (SWE) has been increasingly utilized in medical imaging for its ability to visualize tissue stiffness in vivo, combined with its capability for diagnostic disease assessment [1,2,3,4]. By measuring the speed of the shear waves, SWE provides valuable information about tissue elasticity [5,6] and viscoelasticity [7]. This technology has found applications in the diagnosis of different diseases such as liver fibrosis [8,9,10], breast cancer [11,12], skin lesions [13], and brain diseases such as Alzheimer’s [14].

SWE can be enabled in different modalities, including ultrasound [15,16], magnetic resonance imaging (MRI) [17,18,19], and optical coherence tomography (OCT) [20,21]. These modalities offer complementary advantages in terms of imaging depth, spatial resolution, field of view, and tissue characterization, expanding the capabilities of SWE for clinical applications and research studies.

Ultrasound elastography offers several advantages including almost real-time imaging, low imaging costs, a wide field of view, and significant penetration depth, making ultrasound elastography well-suited for a wide range of clinical applications with good spatial resolution [16,22]. Magnetic resonance elastography (MRE), with its excellent soft-tissue contrast and ability to assess large tissue volumes, provides three-directional and three-dimensional elastograms and is well suited for tissue characterization in various clinical applications, particularly in liver fibrosis assessment and brain studies [23,24]. Optical coherence elastography (OCE), on the other hand, provides high-resolution 3D imaging capabilities and is ideal for applications in ophthalmology and dermatology [25,26].

There are numerous estimators available in elastography; the broad categories of these are outlined in a review paper by Doyley [27]. In shear-wave elastography, there are several approaches including the phase gradient method [28], the time-of-flight method [29], the Fourier estimator [21], the viscoelastic wave inversion method [30], the wave attenuation estimator [31], and different forms of correlation estimators [6,32,33]. In seismography, Aki [32] determined the correlation between transient tremor recordings in one and two dimensions. Later in elastography, Benech and colleagues developed a system of correlating transient shear-wave signals in the framework of time reversal operators [34].

Separately, our reverberant shear-wave elastography employs a narrow-band superposition of random isotropic shear waves propagating in 3D within the tissue, distributed across all directions. A reverberant field is typically produced by multiple shear-wave sources generating waves that reflect and interact with each other, resulting in a complex 3D distribution of shear waves throughout the tissue, with specific limiting ensemble averages described by the real part of the autocorrelation function. Recent studies have demonstrated the effectiveness of reverberant shear-wave elastography in assessing diverse elastic and viscoelastic media [35,36]. However, choosing the appropriate autocorrelation estimator and autocorrelation parameters for a particular elastography application requires careful consideration of various factors, such as the imaging system, tissue properties, experimental conditions, and mechanical excitation [37].

Despite the developments in reverberant SWE, accurately measuring SWS in elastography scans in some situations can be challenging. For example, generating a fully reverberant shear-wave field in 3D typically requires several excitation sources plus reflections, and this is not always achieved in real experiments. Imperfect reverberant shear-wave fields that exhibit a dominant direction or source, especially in proximity to an external source, can be encountered in practice. To address this issue, we introduce a novel autocorrelation estimator based on angular integration for SWS measurement within the general framework established for reverberant shear-wave fields. The effectiveness of the angular integral autocorrelation (AIA) estimator is examined using numerical simulations and experimental data from ultrasound elastography, MRE, and OCE imaging systems.

This paper is organized to first review the fundamental equations for reverberant shear-wave fields and earlier and simple estimators of shear wavelength in these fields. Next, the AIA estimator is introduced. These estimators are then tested in numerical simulations of shear-wave fields and in experiments using ultrasound, MRI, and OCT imaging platforms across a variety of conditions and elastic biomaterials. Final comparisons illustrate the superiority of the proposed AIA estimator, showing its robust performance even in the presence of non-ideal conditions.

2. Theory

When shear waves are generated by multiple excitation sources or wave reflections happen at tissue boundaries, shear-wave interferences generate a reverberant shear-wave field or diffuse shear-wave field [33,38,39]. A fully reverberant shear-wave field can be characterized as the superposition of several planar shear waves propagating in random directions with the same wavenumber $k$ and angular frequency ${\omega }_0$. However, if there is an insufficient number of shear-wave sources or the wave traveling directions are non-random, imperfect, or non-fully reverberant, shear-wave fields are produced. The particle velocity of a fully reverberant shear-wave field $\mathit{V}\left(\mathit{\epsilon },t\right)$ in an isotropic medium can be defined as:

$$\mathit{V}\left(\mathit{\epsilon },t\right)=\sum _{q,l}{\widehat{\mathit{n}}}_{ql}{ v}_{ql} e^{i\left(k{\widehat{\mathit{n}}}_q. \mathit{\epsilon }-{\omega }_{0}t \right)},$$

(1)

where $\mathit{\epsilon }$ represents the position vector, $t$ denotes time, and the indices $q$ and $l$ refer to realizations of the random unit vectors ${\widehat{\mathit{n}}}_q$ and ${\widehat{\mathit{n}}}_{ql}$, respectively. The vector ${\widehat{\mathit{n}}}_q$ represents the random direction of wave propagation, while ${\widehat{\mathit{n}}}_{ql}$ represents a random unit vector indicating the direction of particle motion.${ v}_{ql}$ is an independent, identically distributed random variable representing the magnitude of the particle velocity within a realization. It is worth mentioning that in transversal shear-wave fields, the wave propagation direction is perpendicular to the particle motion, which implies that ${\widehat{\mathit{n}}}_{ql} {. \widehat{\mathit{n}}}_q=0$.

Typically, elastography modalities such as ultrasound elastography, MRE, and OCE are able to measure the particle velocity along a sensor axis (note that MRE can determine displacements along three axes at the expense of longer scan times). It is common practice to designate the sensor axis as the z-axis. Consequently, for the scalar velocity field in the z-direction $V_{\mathit{z}}\left(\mathit{\epsilon },t\right)$ we have:

$$V_{\mathit{z}}\left(\mathit{\epsilon },t\right)=\mathit{V}\left(\mathit{\epsilon },t\right) .{\widehat{\mathit{e}}}_z=\sum _{q,l}{n}_{{ql}_z}{ v}_{ql} e^{i\left(k{\widehat{\mathit{n}}}_q. \mathit{\epsilon }-{\omega }_{0}t \right)},$$

(2)

where ${\widehat{\mathit{e}}}_z$ is the z-direction unit vector and $n_{{ql}_z}={\widehat{\mathit{n}}}_{ql} . {\widehat{\mathit{e}}}_z$ is a scalar random variable. As discussed in the introduction section, for elastography purposes, the SWS or the wave number of the reverberant shear-wave field should be calculated. To calculate the wavenumber $k$ from Equation (2), one approach derives an autocorrelation estimation. So, by considering the autocorrelation function of Equation (2) both in space and time, $B_{V_z{V}_z}\left(∆\mathit{\epsilon }, ∆t\right)$, we derive the following expression:

$$B_{V_z{V}_z}\left(∆\mathit{\epsilon }, ∆t\right)=E\left\{V_{\mathit{z}}\left(\mathit{\epsilon },t\right){V_{\mathit{z}}}^{\mathit{*}}\left(\mathit{\epsilon }+∆\mathit{\epsilon },t+∆t\right)\right\},$$

(3)

where E represents an ensemble average, and the asterisk ($\mathit{\ast }$) denotes the complex conjugate. By substituting Equation (2) into Equation (3), we have:

$$B_{V_z{V}_z}\left(∆\mathit{\epsilon }, ∆t\right)=E\left\{\left(\sum _{q,l}{n}_{{ql}_z}{ v}_{ql} e^{i\left(k{\widehat{\mathit{n}}}_q. \mathit{\epsilon }-{\omega }_{0}t \right)}\right)\times \left(\sum _{q^{\prime },l^{\prime }}{n}_{{q^{\prime }{l}^{\prime }}_z}{ v}_{q^{\prime }{l}^{\prime }} e^{-i\left(k{\widehat{\mathit{n}}}_{q^{\prime }} .\left(\mathit{\epsilon }+∆\mathit{\epsilon }\right)-{\omega }_0\left(t+∆t\right) \right)}\right)\right\}.$$

(4)

Many of the terms in this equation arise from independent realizations and therefore the cross terms are equal to zero. Then, by utilizing the spherical coordinate system described by Aleman-Castañeda et al. [40], Equation (4) simplifies to:

$$B_{V_z{V}_z}\left(∆\mathit{\epsilon },∆t\right)=3{\overline{V_z}}^2 e^{i{\omega }_0∆t}\left\{\frac{{\sin }^2{\theta }_s}{2}\left[j_0\left(k∆\epsilon \right)-\frac{j_1\left(k∆\epsilon \right)}{k∆\epsilon }\right]+\mathrm{c}{\mathrm{os}}^2{\theta }_s \frac{j_1\left(k∆\epsilon \right)}{k∆\epsilon }\right\},$$

(5)

where $\overline{v^2}$ is the ensemble average velocity-squared, ${\theta }_s$ is the angle of $∆\mathit{\epsilon }$ with respect to the sensor axis (considered as the z-axis), $j_0$ and $j_1$ are the first kind spherical Bessel function of zero order and first order, respectively. It is important to note that the vector $∆\mathit{\epsilon }$ is a function of its magnitude $∆\epsilon$ and its angle ${\theta }_s$. The autocorrelation function $B_{V_z{V}_z}\left(∆\mathit{\epsilon },∆t\right)$ clearly depends on the direction of $∆\mathit{\epsilon }$. In simple baseline autocorrelation estimation, $∆\mathit{\epsilon }$ is assumed to align with one of the Cartesian axes; in other words, we have $∆{\epsilon }_x$, $∆{\epsilon }_y$, and $∆{\epsilon }_z$. The angle ${\theta }_s$ for $∆{\epsilon }_x$ and $∆{\epsilon }_y$ is equal to $\pi /2$ and for $∆{\epsilon }_z$ is equal to zero. Therefore, the simple baseline autocorrelation functions are defined as:

$$B_{V_z{V}_z}\left(∆{\epsilon }_x,∆t\right)=\frac{3}{2} {\overline{V_z}}^2 e^{i{\omega }_0∆t} \left[j_0\left(k∆{\epsilon }_x\right)-\frac{j_1\left(k∆{\epsilon }_x\right)}{k∆{\epsilon }_x}\right]$$

(6)

$$B_{V_z{V}_z}\left(∆{\epsilon }_y,∆t\right)=\frac{3}{2} {\overline{V_z}}^2{e}^{i{\omega }_0∆t} \left[j_0\left(k∆{\epsilon }_y\right)-\frac{j_1\left(k∆{\epsilon }_y\right)}{k∆{\epsilon }_y}\right]$$

(7)

$$B_{V_z{V}_z}\left(∆{\epsilon }_z,∆t\right)=3{\overline{V_z}}^2{e}^{i{\omega }_0∆t} \frac{j_1\left(k∆{\epsilon }_z\right)}{k∆{\epsilon }_z}.$$

(8)

The simple baseline autocorrelation estimators are written in terms of spherical Bessel functions, and the wave number $k$ can be calculated using curve fitting. However, it is important to note that when utilizing simple baseline autocorrelation estimators, the autocorrelation data from different directions are ignored, and at minimum, only the autocorrelation function from one direction of lag is considered. This limitation makes the simple baseline autocorrelation estimators highly sensitive to noise and uncertainty, as information from other directions is disregarded. Furthermore, any non-ideal weighting of the shear-wave distribution across all solid angles can cause bias in the estimate from only one or two lag directions of the sampled autocorrelation function.

To address the issue, we introduce a novel approach based on the angular integration of the autocorrelation function. This approach allows us to incorporate the autocorrelation data from various directions while eliminating the dependence of the autocorrelation function on the angle. By performing the angular integral of the autocorrelation function within a two-dimensional plane over the range of 0 to 2π, we obtain the following expression:

$$B_{AI}\left(∆\rho , ∆t\right)=\frac{1}{2\pi }{\int }_0^{2\pi }{B}_{V_z{V}_z}\left(∆\rho , {\theta }_s,∆t\right)d{\theta }_s,$$

(9)

where $B_{AI}\left(∆\rho , ∆t\right)$ represents the AIA function and $∆\rho$ is the one-dimensional shift (or lag) in the autocorrelation argument after integration around ${\theta }_s$. It is important to note that the AIA will reduce to different analytic formulas across different autocorrelation planes. Considering the sensor is aligned with the z-axis, for xz or yz autocorrelation planes, the AIA function can be measured by substituting Equation (5) into Equation (9) and calculating the integral as follows:

$${B_{AI}}_{xz}\left(∆\rho , ∆t\right)={B_{AI}}_{yz}\left(∆\rho , ∆t\right)=\frac{3}{4} {\overline{V_z}}^2 e^{i{\omega }_0∆t} \left[j_0\left(k∆\rho \right)+\frac{j_1\left(k∆\rho \right)}{k∆\rho }\right],$$

(10)

where ${B_{AI}}_{xz}\left(∆\rho , ∆t\right)$ and ${B_{AI}}_{yz}\left(∆\rho , ∆t\right)$ represent AIA functions in the xz and yz planes respectively. These functions account for the autocorrelation data obtained from different directions within each respective plane. However, for the xy autocorrelation plane, the angle ${\theta }_s$ is always equal to $\pi /2$ and so we have the following equation:

$${B_{AI}}_{xy}\left(∆\rho , ∆t\right)=\frac{3}{2} {\overline{V_z}}^2 e^{i{\omega }_0∆t} \left[j_0\left(k∆\rho \right)-\frac{j_1\left(k∆\rho \right)}{k∆\rho }\right],$$

(11)

where ${B_{AI}}_{xy}\left(∆\rho , ∆t\right)$ is the AIA function in the xy plane.

These two descriptions for AIA are independent of the angle ${\theta }_s$, allowing for a straightforward estimation of the wavenumber through curve fitting or other related methods. Moreover, by incorporating all the autocorrelation data from various angles into the calculating autocorrelation function, this method has the potential to be less sensitive to noise and imperfect distributions, resulting in enhanced robustness in SWS measurements.

3. Methods

In this study, we introduce AIA, which calculates the angular integral of the two-dimensional autocorrelations across a designated plane. To assess the efficiency of the AIA approach, some simulations were conducted using the k-Wave simulation toolbox in MATLAB (The MathWorks, Inc. Natick, MA, USA, version 2022b). Subsequently, our proposed method was employed to measure SWS in different modalities, including ultrasound elastography, MRE, and OCE.

In order to provide a comprehensive understanding of the AIA estimator, we detail the algorithm utilized for shear-wave measurements. Figure 1 provides an overview of the algorithm employed to measure the SWS. The case study in this figure is the ultrasound elastography of a breast phantom with a stiff lesion. As depicted in Figure 1, the initial dataset comprises the displacement (or velocity) field obtained from ultrasound elastography, representing the cross-section of an object over time. The displacement field is structured as a 3D real matrix including the depth or axial dimension, lateral dimension, and time dimension. In the first step, an autocorrelation estimation can be applied to the displacement field for each time frame. However, to reduce computational cost and enhance robustness, it is preferable to measure the autocorrelation function in the frequency domain. Therefore, the fast Fourier transform (FFT) is performed in the time domain, followed by the selection of the maximum intensity in the frequency domain, corresponding to the excitation frequency. The resulting data in this step are a 2D complex matrix, comprising both magnitude and phase components, as presented in Figure 1.

After applying filtering techniques including finite impulse response (FIR) filtering, 2D bandpass filtering, and median filtering to experimental data to reduce unwanted noise in the results and enhance the quality and reliability of the data, a 2D autocorrelation window with a defined window size and window step (jump) moves through the 2D matrix. The autocorrelation of each window is computed, and the real part of the results is considered for further analysis. As illustrated in Figure 1, each autocorrelation window generates a 2D real matrix, where the $∆\mathit{\epsilon }$ values within this matrix are in Cartesian coordinates but can be converted to polar $(r$, $\theta )$ values using ${\epsilon }_x=r\cos \theta$ and ${\epsilon }_y=r\sin \theta$. To calculate the angular integral, we sum over the autocorrelation matrix across constant values of $r$, resulting in a one-dimensional autocorrelation function of $∆\rho$. Full curve fitting or other related methods can be utilized to estimate the unknown wavenumber k or the SWS for each measured curve. In this study, MATLAB curve fitting was utilized. This approach enables an increase in the accuracy and reliability of the estimates for the SWS as it considers the autocorrelation data from all angles and reduces the impact of noise and uncertainty. Finally, the SWS of the entire 2D cross-section is determined by combining the estimated SWS from all autocorrelation windows across the entire cross-section.

This algorithm can also be extended for SWS estimation in a 3D medium by dividing the 3D medium into 2D planes. The same approach can be applied to each plane, treating them as individual cross-sections, and then the SWS can be combined to obtain a 3D representation of the SWS distribution, including possible anisotropy [40]. This adaptation allows for the assessment of SWS variations in different planes of the 3D medium, providing a more comprehensive understanding of the mechanical properties throughout the entire volume. In this study, alongside the comparison to ground truth, the SNR within homogeneous regions is calculated using the following equation [28].

$$\mathrm{S}\mathrm{N}\mathrm{R}=10 {\log }_{10} \left(\frac{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e} \mathrm{o}\mathrm{f} \mathrm{s}\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{r}-\mathrm{w}\mathrm{a}\mathrm{v}\mathrm{e} \mathrm{s}\mathrm{p}\mathrm{e}\mathrm{e}\mathrm{d}}{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{d} \mathrm{d}\mathrm{e}\mathrm{v}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \mathrm{o}\mathrm{f} \mathrm{s}\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{r}-\mathrm{w}\mathrm{a}\mathrm{v}\mathrm{e} \mathrm{s}\mathrm{p}\mathrm{e}\mathrm{e}\mathrm{d}} \right)$$

(12)

This metric can be used across each of the examples and modalities presented in the results and is computed using the largest region of interest (ROI) that can be selected within each dataset within the nominally homogeneous background or interior. The SNR measurement for background regions (homogeneous regions) is chosen because it could be applied to all our examples from all modalities, providing a consistent measure across all examples. Ideally, it should achieve high levels, characterized by a uniform mean and low standard deviation across a homogeneous background in low noise.

3.1. Numerical Simulations

To demonstrate the effectiveness of the AIA estimator for SWS estimation, a simulation was conducted in MATLAB (http://www.k-wave.org/documentation.php, accessed on 1 January 2024). The simulation aimed to model real experimental conditions for different elastography modalities by considering various aspects of wave generation, shear-wave propagation, and shear-wave interactions in a medium.

As illustrated in Figure 2a, the simulated model consisted of a stiff branching y-shaped tube placed within a uniform background material. The model was represented as a cube with dimensions of 120 mm × 120 mm × 120 mm, and the background material was assumed to be a uniform isotropic material with a SWS of 1 m/s. The y-shaped tube was also assumed to be a uniform isotropic material with a cross-section radius of 15 mm and a SWS of 3 m/s. The density of the medium, both for the background and the y-shaped tube, was specified as 1000 $\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$. In the simulation, the viscoelastic properties of the materials were defined using the classical Kelvin–Voigt absorption model. This model considers the absorption coefficient to be proportional to the square of the frequency in the low-frequency limit [41,42]. For both the background material and the y-shaped tube, the absorption coefficient was set to 0.5 $\mathrm{d}\mathrm{B}/\left({\mathrm{M}\mathrm{H}\mathrm{z}}^2\mathrm{c}\mathrm{m}\right)$. By incorporating the viscoelastic properties and absorption coefficients based on the Kelvin–Voigt model, the simulation is designed to capture the realistic behavior of shear-wave propagation and absorption in the simulated materials.

The shear waves were generated by the k-Wave point velocity excitation sources placed in random positions near the boundaries of the simulated medium. This approach ensured the generation of a fully reverberant shear-wave field. To verify that the excitation sources did not directly generate waves within the y-shaped tube, the length of the tube was made smaller than the overall dimensions of the cube. This arrangement allowed the waves to be generated solely within the cube and then propagate both inside the cube and the tube. This setup reflects the typical scenario in many elastography experiments, where the waves are generated outside the region of interest (ROI) and propagate into the target area.

A white Gaussian noise with a signal-to-noise ratio (SNR) of 10 dB was introduced into the shear-wave field in the k-Wave simulation. The addition of this noise in the simulation enables the evaluation of the robustness and performance of the AIA approach under more realistic conditions. The impact of noise on the accuracy performance of AIA is further examined in Appendix A.

Furthermore, it was assumed that in this elastography experiment, only the velocity field along the z-axis could be measured. This assumption aligns with imaging modalities where the sensitive displacement measurements are mainly directed along one specific direction. By including these considerations in the simulation, a realistic experimental setup that closely resembles the conditions encountered in practical elastography studies was created. The particle velocity field along the z-axis, after sufficient time (0.133 s in this case) for steady-state wave interaction and the generation of a reverberant shear-wave field, is shown in Figure 2b. Figure 2c presents a focused view of the particle velocity field within the y-shaped tube. It is clear in this figure that a fully reverberant shear-wave field is generated in both the background and the y-shaped tube. The distinct difference between the wave amplitudes and wavelengths in a y-shaped tube and the background makes the y-shaped tube recognizable in the shear-wave field. To reduce computational costs, the 3D ROI has been selectively cut.

3.2. Ultrasound Elastography

In this section, the effectiveness of the AIA estimator for estimating the SWS in ultrasound elastography is investigated. To provide a comprehensive assessment across various experimental scenarios and conditions, two different ultrasound experiments with different excitation systems and field-of-view configurations were studied. In the first study, a series of ultrasound elastography experiments were performed utilizing an elastic breast phantom containing a stiff lesion. These experiments were conducted at different excitation frequencies. In these experiments, imperfect reverberant shear-wave fields were generated using two excitation sources positioned on opposite sides of the phantom. The second study focused on the utilization of datasets obtained from in vivo ultrasound elastography of the human liver–kidney region with a fully reverberant shear-wave field. In this case, data were collected at different excitation frequencies.

3.2.1. Breast Phantom Ultrasound Elastography

In the first ultrasound elastography study, a CIRS breast phantom (model 509, CIRS Inc., Norfolk, Virginia, USA) was utilized. The CIRS breast phantom is designed to simulate human breast tissue characteristics and closely replicates the mechanical properties of actual breast tissue. The breast phantom contains several lesions of different sizes. In this study, to illustrate the capability of the AIA estimator to differentiate the stiffer regions from the background material (the region of interest), the focus was directed toward a specific 10 mm diameter lesion within the phantom.

To generate a reverberant shear-wave field, two mechanical vibration sources (model 4810, Bruel and Kjaer, Naerum, Denmark) and two miniature vibration sources (model NCM02-05-005-4 JB, H2W, Linear Actuator, Santa Clara, CA, USA) were employed in contact with the breast phantom. These vibrators were positioned at four points, forming a rectangle around the phantom’s sides. The excitation signal for the two mechanical vibration sources was supplied by a power amplifier (model 2718, Bruel and Kjaer, Naerum, Denmark), while for the two miniature vibration sources, it was provided by a digital power amplifier (model LP-2020 A+, Lepai, Bukang, China) driven by a dual-channel function generator (model AFG3022B, Tektronix, Beaverton, OR, USA). The ultrasound elastography experiment was performed at various excitation frequencies, including 900 Hz, 600 Hz, and 400 Hz.

The ultrasound elastography experiment utilized a Verasonics ultrasound system (V-1, Verasonics, Kirkland, WA, USA). This system enables high frame-rate acquisition and a coherent plane-wave compounding acquisition scheme. The Verasonics system was connected to a linear-array ultrasound transducer (Model L7-4, ATL, Bothell, WA, USA). To track the induced displacements caused by shear-wave propagation, a Loupas estimator was employed [43]. In the ultrasound elastography experiment, a 3D matrix of in-phase and quadrature (IQ) data was collected and stored for subsequent postprocessing. The IQ data contain information about the amplitude and phase of the ultrasound signals received by the transducer. To calculate the axial particle displacements, a frame-to-frame analysis was performed on the acquired 3D IQ data. Then, the amplitude and phase of displacement at each pixel were estimated. In the ultrasound elastography experiment, a 3D matrix of in-phase and quadrature (IQ) data was collected and stored for subsequent postprocessing. The IQ data contain information about the amplitude and phase of the ultrasound signals received by the transducer. To calculate the axial particle displacements, frame-to-frame analysis was performed on the acquired 3D IQ data. Then, the amplitude and phase of displacement at each pixel were estimated. In all the experiments, the center frequency of the ultrasound transducer was set to 5 MHz. The tracking pulse repetition frequency (PRF) was adjusted to acquire at least 20 samples per cycle of the vibration frequency. Further information regarding the characteristics of the phantom, the parameters of the ultrasound elastography experiment, and data postprocessing to determine the displacement field can be found in Ormachea and Parker [28,33].

Figure 3a displays the B-mode ultrasound scan of the breast phantom with a 10 mm diameter lesion positioned at the center. In this experimental data, shear waves were generated within the phantom at various excitation frequencies, and ultrasound elastography scans were performed during the excitation process. The phase map of the shear-wave field at the frequency of 900 Hz is presented in Figure 3b. This representation describes an imperfect reverberant shear-wave field within the phantom that is more bi-directional. As evident in the figure, the lesion is visible in the B-mode scan, while it cannot be recognized in the shear-wave field.

3.2.2. Liver–Kidney Ultrasound Elastography

In the next study, to assess the effectiveness of the AIA estimator, datasets derived from ultrasound elastography of the human liver–kidney region were investigated. The liver scan was obtained in conjunction with the small study reported by Ormachea et al. [44], which was performed under the requirements of informed consent of the Southwoods Imaging Clinical Institutional Review Board. A Verasonics ultrasound system (Vantage-128™, Verasonics, Kirkland, WA, USA), connected to a convex ultrasound probe (model C4-2, ATL, Bothell, WA, USA), was used with a 3 MHz center frequency. To induce the reverberant shear-wave field, a custom-made portable trifold futon measuring 70 cm × 60 cm × 10 cm was utilized. This setup incorporated eight electromagnetic vibration drivers (Quad Resonator Model EI718™, Elastance Imaging LLC, Columbus, OH, USA) positioned within a square configuration. These drivers were designed to transmit displacements and shear waves into the body surrounding the liver and kidney. Each driver, measuring 4 inches in diameter, could be powered by standard audio amplifiers. These datasets include the B-mode ultrasound scan and displacement field data from fully reverberant shear waves at a broad range of excitation frequencies, including 702 Hz, 585 Hz, 468 Hz, 351 Hz, 234 Hz, and 117 Hz. Additional details about the data postprocessing to obtain the displacement field are provided in Ormachea et al. [44].

Figure 4a displays the B-mode ultrasound scan, which provides an anatomical view of the liver and kidney. In the scan, the different tissue layers can be observed, including the abdominal muscle and fat layers (zone A), the distinct darker zones representing the liver (zone L), and the kidney (zone K). Figure 4b presents the phase map of the shear-wave field at the frequency of 702 Hz. This representation depicts a randomized, reverberant shear-wave field throughout the abdominal region.

3.3. Magnetic Resonance Elastography

In the next step, an MRE experiment was conducted on a brain-mimicking phantom with two stiff lesions. Figure 5a displays the brain-shaped phantom, along with the coordinate system for reference. The brain-shaped phantom was constructed using a suspension composed of bovine gelatin (200 bloom: Sigma Aldrich Chemicals, St. Louis, MO, USA), de-ionized water (18 MΩ), and ethylenediamine tetra-acetic acid (Sigma Aldrich Chemicals, St. Louis, MO, USA). A controlled and reproducible process was used (as described by Doyley et al. [45]) to fabricate the phantom with precise dimensions of 180 mm in the long axis, 130 mm in the short axis, and 70 mm in height. The brain phantom contained two spherical gelatin lesions, measuring 18 mm and 12 mm in diameter, respectively. The gelatin concentration in the background region was 8%, while both lesions had a gelatin concentration of 18%. The remaining composition for the background consisted of 92% water, while the lesion composition consisted of 81.64% water and 0.36% copper sulfate. The addition of copper sulfate served to provide contrast between the lesion and the background, enhancing the anatomical image in the scan. The fabrication process ensured the creation of a brain-shaped phantom with accurately defined dimensions and distinct gelatin concentrations, allowing for reliable experimental testing and imaging analysis.

The single-shot echo-planar imaging (EPI) sequence was used to compute the time-varying harmonic tissue displacements. The MR scanner was configured with an echo time of 76.0 ms and a repetition time of 8640 ms. A total of forty axial slices with an isotropic voxel size of 1.6 mm were obtained for the phantom, covering a field-of-view of 153 mm × 153 mm × 60 mm. The shear-wave field was generated within the phantom at the frequency of 200 Hz by a pneumatic mechanical actuator with a passive driver (Resoundant, Inc., Rochester, MN, USA). The mechanical vibrator was positioned at the bottom of the brain phantom and operated in the y-direction. All elastographic imaging procedures were conducted using a whole-body 3T MRI scanner (Prisma, Siemens, Erlangen, Germany) equipped with a 20-channel head coil. Each 3D dataset acquisition required approximately six minutes. Data collection in MRE was carried out with a 1.6 mm isotropic voxel size. The single-shot echo-planar imaging (EPI) sequence [46] was utilized to measure the resulting time-varying harmonic tissue displacements. More details about the phantom properties, MRE experiment parameters, and data postprocessing to determine displacement field are described in Kabir et al. [47].

In the MRI scan of the brain phantom (see Figure 5b), two spherical lesions are clearly visible. The darker small regions visible in the MRI scan are air bubbles. The phase map of X-motion and Z-motion displacement fields at the frequency of 200 Hz are displayed in Figure 5c,d, respectively. It is worth noting that X-motion and Z-motion displacement fields were exclusively employed as they exhibited the strongest signal. As shown in Figure 5, shear-wave fields with the dominant directions are generated in both X-motion and Z-motion. The presence of lesions is visible within these displacement fields as longer spatial wavelengths and the bending of wavefronts around the lesions.

3.4. Optical Coherence Elastography

To assess the effectiveness of the AIA approach in OCE, an OCE experiment was conducted on a gelatin phantom featuring a cylindrical stiff lesion. For this experiment, a custom-built OCT system equipped with a swept-source laser with a center wavelength of 1310 nm and a bandwidth of 140 nm (HSL-2100-HW, Santec, Aichi, Japan) was employed. In the air, the system achieved maximum lateral and depth resolutions of 20 μm and 6 μm, respectively, with a sensitivity of approximately 110 dB. For this experiment, the field of view was set at 10 mm × 10 mm, with 100 scanning points in each direction on the phantom. Consequently, the lateral scanning resolution of the OCE experiment was adjusted to 100 μm × 100 μm. The penetration depth of the OCE experiment was measured to be 0.5 mm. Synchronization between the swept-source optical coherence tomography (SS-OCT) and a mechanical excitation system was achieved using LabVIEW (Version 14, National Instruments, Austin, TX, USA).

The phantom and its lesion were composed of 5% gelatin concentration and 10% gelatin powder concentration, respectively, designed to exhibit distinct mechanical properties and SWS characteristics. To introduce the light scattering properties to the phantom, a 3% concentration of intralipid powder (coffee creamer) was added to both the background and lesion. The remaining composition for the background consisted of 1% salt and 91% water, while the lesion composition comprised 1% salt and 86% water.

A reverberant shear-wave field at the frequency of 1500 Hz was induced within the phantom using a custom 3D-printed multi-pronged ring connected to a piezoelectric actuator (BA4510, PiezoDrive, Callaghan, NSW, Australia). The custom ring featured eight concentric points equidistant from the center, with a 10 mm diameter aperture to accommodate SS-OCT scans. The custom ring was placed gently onto the gelatin phantom to induce reverberant shear waves. Excitation signals were provided using a function generator (AFG320, Tektronix, Beaverton, OR, USA) and an amplifier (PDu150, PiezoDrive, Callaghan, New South Wales, Australia). Additional information regarding the optical setup and excitation system used in OCE can be found in Ge et al. [48]. The acquisition approach and data postprocessing developed by Zvietcovich et al. [6] were employed to obtain the displacement field.

A three-dimensional OCT scan along with the corresponding displacement field of shear waves for the phantom containing a stiff lesion at the frequency of 1500 Hz is presented in Figure 6. Notably, a pie-cut is applied to view the interior. In the three-dimensional OCT scan (see Figure 6a), the presence of the cylindrical lesion is visually revealed through the height difference between the phantom and the lesion at the center. However, it should be emphasized that both the phantom and its lesion were intentionally designed with identical scattering properties. Therefore, despite differences in surface height, the 3D OCT scan does not distinguish between the stiff lesion and background due to their matching scattering properties. Figure 6b clearly illustrates the generation of circular waves at the boundaries of the phantom, which subsequently propagate and interfere with each other within the internal region. This results in a wave field that exhibits a more directional nature at the outer phantom boundaries and a more reverberant behavior within the internal region.

4. Results and Discussion

4.1. Numerical Simulation

In the k-Wave simulation, access to the complete 3D dataset allowed for the analysis of shear-wave behavior in various planes, including the xy plane, xz plane, and yz plane. Within each plane, measurements were conducted for the simple baseline autocorrelation functions in two directions and the AIA function. Specifically, for the xy plane, measurements involved the simple baseline autocorrelation function along the x-axis, denoted as $B_{V_z{V}_z}\left({∆\epsilon }_x,∆t\right),$ and the y-axis, denoted as $B_{V_z{V}_z}\left({∆\epsilon }_y,∆t\right)$, as well as the AIA in the xy plane, denoted as ${B_{AI}}_{xy}\left(∆\rho , ∆t\right)$.

Since we assumed that the detected velocity field is aligned with the z-axis, the presumed imaging direction and the autocorrelation function in the xy plane remained consistent regardless of the autocorrelation direction. In other words, the simple baseline autocorrelation functions yield similar results for both $B_{V_z{V}_z}\left({∆\epsilon }_x,∆t\right)$ and $B_{V_z{V}_z}\left({∆\epsilon }_y,∆t\right)$ directions, and ${\theta }_s$ in Equation (5) is equal to $\pi /2$ for both ${∆\epsilon }_x$ and ${∆\epsilon }_y$, resulting in Equations (6) and (7). Similarly, measurements of the simple baseline autocorrelation function were conducted along the x-axis $B_{V_z{V}_z}\left({∆\epsilon }_x,∆t\right)$ and the z-axis $B_{V_z{V}_z}\left({∆\epsilon }_z,∆t\right)$ in the xz plane as well as along the y-axis $B_{V_z{V}_z}\left({∆\epsilon }_y,∆t\right)$ and the z-axis $B_{V_z{V}_z}\left({∆\epsilon }_z,∆t\right)$ in the yz plane. In both the xz and yz planes, the angle ${\theta }_s$ for ${∆\epsilon }_x$ and ${∆\epsilon }_y$ is equal $\pi /2,$ and the simple baseline autocorrelation functions follow Equations (6) and (7). For ${∆\epsilon }_z$, ${\theta }_s$ is zero and the simple autocorrelation function $B_{V_z{V}_z}\left({∆\epsilon }_z,∆t\right)$ follows Equation (8). It is important to note that when measuring simple baseline autocorrelation functions in other directions with different ${\theta }_s$, Equation (5) should be employed.

Through the AIA estimation conducted in each xy, xz, and yz plane (per Equations (10) and (11)), a comprehensive analysis of the shear waves in different planes was conducted, providing valuable insights for SWS estimation and tissue characterization. Figure 7 presents the SWS measured from the k-Wave simulation at a frequency of 200 Hz, using the simple baseline autocorrelation approaches in different planes and directions (left column) and our proposed AIA approach in different planes (middle column). The autocorrelation window size in all measurements was 15 mm × 15 mm, and the window step was 3 mm.

Table 1. SWS in the background material and the branching tube estimated using simple autocorrelation approaches and the AIA approach in k-Wave simulation across different autocorrelation planes.

Autocorrelation Plane	SWS (m/s) in Background Using Simple Autocorrelation (First Axis)	SWS (m/s) in Background Using Simple Autocorrelation (Second Axis)	SWS (m/s) in Background Using AIA	SWS (m/s) in Branching Tube Using Simple Autocorrelation (First Axis)	SWS (m/s) in Branching Tube Using Simple Autocorrelation (Second Axis)	SWS (m/s) in Branching Tube Using AIA
xy plane	1.12 ± 0.14	1.03 ± 0.09	1.02 ± 0.04	2.90 ± 0.74	2.70 ± 0.42	2.98 ± 0.57
xz plane	1.08 ± 0.23	1.08 ± 0.13	1.14 ± 0.12	2.93 ± 0.61	2.51 ± 0.47	3.19 ± 0.74
yz plane	1.07 ± 0.17	0.97 ± 0.08	1.05 ± 0.08	2.81 ± 0.68	2.50 ± 0.36	3.08 ± 0.65

showcases a comparison of the estimated SWS using simple autocorrelation approaches and the AIA approach in the background material and the branching tube across different autocorrelation planes. It is worth noting that the SWS and SNR for all measurements were computed using the almost-full 3D background and the y-shaped tube, with some distance to the tube-background boundaries. The average error for estimating background SWS was 7% using either simple autocorrelation approaches or AIA. However, the average error for estimating SWS in the branching tube using simple autocorrelation approaches was 9%, and using the AIA approach, the average error was 3%. Another important parameter in elastography estimation is the ratio of tube SWS to background SWS.

Table 2. SWS ratio between the branching tube and the background material estimated using simple autocorrelation approaches and the AIA approach in k-Wave simulation across different autocorrelation planes.

Autocorrelation Plane	SWS Ratio Using Simple Autocorrelation (First Axis)	SWS Ratio Using Simple Autocorrelation (Second Axis)	SWS Ratio Using AIA	SWS Ratio Estimation Error Using Simple Autocorrelation (First Axis)	SWS Ratio Estimation Error Using Simple Autocorrelation (Second Axis)	SWS Ratio Estimation Error Using AIA
xy plane	2.59 ± 0.74	2.62 ± 0.47	2.92 ± 0.57	14%	13%	3%
xz plane	2.71 ± 0.81	2.32 ± 0.52	2.80 ± 0.71	10%	23%	7%
yz plane	2.63 ± 0.76	2.58 ± 0.43	2.93 ± 0.66	12%	14%	2%

presents a comparison between the estimated SWS ratio using simple autocorrelation approaches and the AIA approach in different autocorrelation planes. The average error of the estimated SWS ratio using simple autocorrelation approaches is 14%, whereas using the AIA approach reduces this error to 4%.

Table 3. SNR of estimated SWS in the background material for k-Wave simulation using simple autocorrelation approaches and the AIA approach across different autocorrelation planes.

Autocorrelation Plane	SNR (dB) in Background Using Simple Autocorrelation (First Axis)	SNR (dB) in Background Using Simple Autocorrelation (Second Axis)	SNR (dB) in Background Using AIA	AIA Improvement over Simple Autocorrelation (First Axis)	AIA Improvement over Simple Autocorrelation (Second Axis)
xy plane	8.94	10.61	13.69	53%	29%
xz plane	6.78	9.13	9.68	43%	6%
yz plane	8.12	10.69	10.98	35%	3%

presents the SNR of estimated SWS in the background material for k-Wave simulation using simple autocorrelation approaches and the AIA approach along with the improvement achieved by the AIA approach in different autocorrelation planes (xy, xz, and yz planes).

The results clearly demonstrate that the AIA approach significantly enhances the SNR in estimating SWS in the background, with an average improvement of 28% at the frequency of 200 Hz. Moreover, we can further improve the estimation by calculating the median, or averaging the estimated SWS in the xy, xz, and yz autocorrelation planes, as shown in the right column of Figure 7. The application of averaging leads to an SNR increase of more than 43%, while employing the median increases the SNR by more than 45% at the excitation frequency of 200 Hz. These findings highlight the effectiveness of the AIA approach and its potential to enhance the accuracy and reliability of SWS estimation.

Extensive comparisons of the SWS ratio and SNR using the average of simple autocorrelation approaches and the AIA approach across different autocorrelation planes are explored in

Table A1. Comparison of SWS in the background material and the branching tube, and the SWS ratio estimated using the average of simple autocorrelation approaches across different autocorrelation planes.

Autocorrelation Plane	SWS (m/s) in Background Using an Average of Simple Autocorrelation Approaches	SWS (m/s) in Branching Tube Using an Average of Simple Autocorrelation Approaches	SWS Ratio Using an Average of Simple Autocorrelation Approaches	SWS Ratio Estimation Error Using an Average of Simple Autocorrelation Approaches	SWS Ratio Estimation Error Using AIA
xy plane	1.08 ± 0.08	2.80 ± 0.41	2.59 ± 0.43	14%	3%
xz plane	1.08 ± 0.12	2.72 ± 0.41	2.52 ± 0.47	16%	7%
yz plane	1.02 ± 0.08	2.66 ± 0.35	2.61 ± 0.40	13%	2%

and

Table A2. SNR comparison of estimated SWS in the background material using the average of simple autocorrelation approaches and the AIA approach across different autocorrelation planes.

Autocorrelation Plane	SNR (dB) in Background Using the Average of Simple Autocorrelation Approaches	SNR (dB) in Background Using AIA	AIA Improvement over Simple Autocorrelation
xy plane	11.40	13.69	20%
xz plane	9.53	9.68	2%
yz plane	10.80	10.98	2%

in Appendix A. Furthermore, the effect of different levels of added noise on the SWS estimated using the AIA approach across different autocorrelation planes is investigated through Figure A2 and Figure A3 and

Table A3. Comparison of SWS both in the background material and the branching tube estimated using the AIA approach in different levels of added Gaussian white noise and different autocorrelation planes.

Added Gaussian White Noise	Autocorrelation Plane	SWS (m/s) in Background Using AIA	SWS (m/s) in Branching Tube Using AIA	SWS Ratio Using AIA	SWS Ratio Estimation Error Using AIA	SWS Ratio Estimation Error Using an Average of Simple Autocorrelation Approaches
	xy plane	1.03 ± 0.05	2.92 ± 0.55	2.83 ± 0.55	6%	17%
1 dB	xz plane	1.15 ± 0.12	3.01 ± 0.72	2.62 ± 0.68	13%	21%
	yz plane	1.06 ± 0.08	2.95 ± 0.57	2.78 ± 0.58	7%	17%
	xy plane	1.03 ± 0.06	2.89 ± 0.50	2.81 ± 0.51	6%	18%
−1 dB	xz plane	1.15 ± 0.12	3.09 ± 0.73	2.69 ± 0.69	10%	22%
	yz plane	1.06 ± 0.08	2.90 ± 0.53	2.74 ± 0.54	9%	17%
	xy plane	1.03 ± 0.07	2.38 ± 0.59	2.31 ± 0.59	23%	37%
−10 dB	xz plane	1.13 ± 0.11	1.98 ± 0.88	1.75 ± 0.80	42%	52%
	yz plane	1.05 ± 0.07	2.16 ± 0.71	2.06 ± 0.69	31%	41%

in Appendix A.

4.2. Ultrasound Elastography

In ultrasound elastography, the dataset primarily comprises the z-oriented displacement field within the imaged xz plane. These data reflect the propagation and interaction of shear waves within the examined tissue. To estimate the SWS map, a simple autocorrelation approach in the x-direction, a simple autocorrelation approach in the z-direction, and the AIA approach were employed.

4.2.1. Breast Phantom Ultrasound Elastography

In our study, different autocorrelation approaches were applied to estimate SWS in ultrasound elastography of the breast phantom with a lesion, conducted under conditions of imperfect reverberant shear-wave fields at various excitation frequencies, including 400 Hz, 600 Hz, and 900 Hz. For these measurements, a square autocorrelation window of dimensions 7.5 mm and a window step of 0.5 mm were utilized.

Figure 8 displays the estimated SWS using simple baseline autocorrelation approaches in the x-direction (see Figure 8a) and z-direction (see Figure 8b), as well as the AIA approach (see Figure 8c) for the excitation frequency of 900 Hz. As shown in Figure 3b, the shear-wave field exhibits an imperfect reverberant nature with prominently directional waves in the x-direction. Under these circumstances, the baseline autocorrelation in the z-direction measures a longer wavelength, resulting in an overestimation of SWS, as evident in Figure 8b. Consequently, the baseline autocorrelation in the z-direction does not yield satisfactory results for SWS estimation; even the difference in SWS between the phantom and lesion is undetectable, as seen in Figure 8b.

The mean SWS and the standard deviation for the CIRS breast phantom were reported as 2.27 ± 0.19 m/s in Ormachea and Parker [28].

Table 4. SWS in the background material for ultrasound elastography of the breast phantom estimated using simple autocorrelation approaches (kx and kz) and the AIA approach across different excitation frequencies, errors were calculated using the SWS value of 2.27 (m/s) for the breast phantom as obtained in Ormachea and Parker [28].

Excitation Frequency	SWS (m/s) in Background Using Simple Autocorrelation (kx)	SWS (m/s) in Background Using Simple Autocorrelation (kz)	SWS (m/s) in Background Using AIA	SWS Estimation Error Using Simple Autocorrelation (kx)	SWS Estimation Error Using Simple Autocorrelation (kz)	SWS Estimation Error Using AIA
900 Hz	2.49 ± 0.17	3.74 ± 1.09	2.45 ± 0.11	10%	65%	8%
600 Hz	2.34 ± 0.13	3.36 ± 0.91	2.22 ± 0.11	3%	48%	2%
400 Hz	2.17 ± 0.13	3.32 ± 0.75	2.10 ± 0.12	4%	46%	7%

presents the estimated SWS in the background using simple autocorrelation approaches (kx and kz) and the AIA approach across different excitation frequencies. The average SWS and SNR for simple autocorrelation approaches and AIA were calculated using the almost-full 2D background, with some distance to the lesion boundaries. The average SWS estimated using simple autocorrelation kx (parallel to the wave direction) and AIA at different excitation frequencies closely align with the value reported by Ormachea and Parker [28]. However, as anticipated, the simple autocorrelation kz (perpendicular to the wave direction) estimates notably higher SWSs for the background at different excitation frequencies. Considering the reported value of 2.27 ± 0.19 m/s as the ground truth, the SWS estimation error for simple autocorrelation approaches and AIA are calculated and presented in Table 4. The average error for SWS estimation using AIA is 6%, while it is 6% and 53% using kx and kz, respectively.

Table 5. SNR of estimated SWS in the background material for the breast phantom ultrasound elastography using simple autocorrelation approaches and the AIA approach across different excitation frequencies.

Excitation Frequency	SNR (dB) in Background Using Simple Autocorrelation (kx)	SNR (dB) in Background Using Simple Autocorrelation (kz)	SNR (dB) in Background Using AIA	AIA Improvement over kx	AIA Improvement over kz
900 Hz	11.73	5.36	13.38	14%	150%
600 Hz	12.49	5.66	13.01	4%	130%
400 Hz	12.35	6.43	12.44	1%	93%

presents the SNR of estimated SWS within the background material for ultrasound elastography of the breast phantom, employing simple autocorrelation approaches (kx and kz) as well as the AIA approach at different excitation frequencies. The estimated SWS using AIA demonstrates a major 150% improvement in SNR within the background region when compared to the kz approach. The SWS obtained using the AIA approach also displays an enhancement in estimation accuracy compared to the simple autocorrelation approach in the x-direction (kx). This improvement is evidenced by a 14% increase in the SNR within the background region. Similar outcomes were obtained at excitation frequencies of 600 Hz and 400 Hz. The presented results in Figure 8 and both Table 4 and Table 5 indicate that AIA is effective in estimating SWS, even in the presence of a highly directional wave field.

4.2.2. Liver–Kidney Ultrasound Elastography

Figure 9a,b presents the SWS obtained from simple baseline autocorrelation functions in the x-direction and z-direction at the frequency of 702 Hz, while Figure 9c illustrates the SWS obtained using the AIA approach at the frequency of 702 Hz. All measurements depicted in Figure 9 were conducted using a square autocorrelation window of dimension 18.5 mm and a window step of 1.1 mm. To evaluate the accuracy of the AIA approach, the SWS elastogram at the frequency of 702 Hz is associated with the ultrasound B-mode scan in Figure 9d. As evident, the estimated SWS using the AIA estimator effectively outlined the structural details of tissues within the liver–kidney region including the abdominal muscle and fat layers (zone A), the liver (zone L), and the kidney (zone K), resembling the structures observed in the B-mode scan.

Consistent outcomes were achieved across a range of excitation frequencies, including 585 Hz, 468 Hz, 351 Hz, 234 Hz, and 117 Hz. In order to analyze the data and compute the SWS and SNR, the largest region of interest (ROI) that can be selected within the liver and kidney were considered for simple autocorrelation approaches and AIA. The frequency-resolved mean and standard deviation values for the liver and the kidney are presented in

Table 6. SWS in the liver and the kidney estimated using simple autocorrelation approaches and the AIA approach in liver–kidney ultrasound elastography across different excitation frequencies.

Excitation Frequency	SWS (m/s) in Liver Using Simple Autocorrelation (kx)	SWS (m/s) in Liver Using Simple Autocorrelation (kz)	SWS (m/s) in Liver Using AIA	SWS (m/s) in Kidney Using Simple Autocorrelation (kx)	SWS (m/s) in Kidney Using Simple Autocorrelation (kz)	SWS (m/s) in Kidney Using AIA
702 Hz	3.13 ± 0.12	4.07 ± 0.12	2.69 ± 0.08	4.93 ± 0.42	4.47 ± 0.25	3.98 ± 0.22
585 Hz	3.12 ± 0.54	3.12 ± 0.15	2.34 ± 0.25	3.81 ± 0.37	3.30 ± 0.37	2.98 ± 0.14
468 Hz	2.71 ± 0.27	2.72 ± 0.11	1.92 ± 0.20	3.11 ± 0.18	2.56 ± 0.10	2.48 ± 0.10
351 Hz	2.47 ± 0.20	2.05 ± 0.12	1.96 ± 0.11	3.25 ± 0.90	2.42 ± 0.27	2.60 ± 0.28
234 Hz	1.61 ± 0.14	2.02 ± 0.10	1.60 ± 0.09	2.18 ± 0.39	1.37 ± 0.05	1.66 ± 0.09
117 Hz	0.80 ± 0.28	0.80 ± 0.05	0.71 ± 0.08	0.95 ± 0.29	0.89 ± 0.05	0.88 ± 0.08

.

Figure 10a,b represents the SNR of measured SWS, employing simple autocorrelation approaches (kx and kz) and the AIA approach within the liver and kidney regions, respectively. The results illustrate that the SNR for simple autocorrelation in the z-direction (kz) is generally higher than the SNR for AIA. However, an examination of the estimated SWS ratio between the liver and kidney in Figure 10c reveals that the SWS ratio estimated by kz consistently remains around 1.0 at different excitation frequencies. This indicates that kz consistently estimates nearly the same SWS for both the liver and kidney, which is inaccurate, as the kidney is commonly found to be stiffer than the liver [49,50]. Furthermore, as evident in Figure 9b, the estimated SWS using kz fails to visualize the structural details of tissues within the liver–kidney region. These results emphasize the effectiveness of the AIA approach in accurately estimating and visualizing the SWS in ultrasound elastography of the in vivo liver–kidney region.

4.3. Magnetic Resonance Elastography

To assess the capability of AIA in visualizing the lesions in MRE scans, our proposed estimator was applied to the 3D MRE datasets obtained from the brain-mimicking phantom experiment. By leveraging the autocorrelation approaches, our objective was to extract the SWS map and use it for lesion visualization and characterization. This study specifically focused on the X-motion and Z-motion signals, which exhibited the strongest signal strengths in this configuration.

Utilizing these MRE experimental data, SWS maps were estimated by employing simple autocorrelation approaches in the x and y directions, as well as the AIA approach. Figure 11 presents the estimated SWS within the brain-mimicking phantom at an excitation frequency of 200 Hz. Figure 11a,b displays the estimated SWS for the X-motion MRE data using the simple baseline autocorrelation in the x-direction and y-direction, respectively. Figure 11c showcases the estimated SWS using AIA for X-motion MRE data. Similarly, Figure 11d–f represents the estimated SWS for the Z-motion MRE data using the simple baseline autocorrelation in the x and y directions, and the AIA approach, respectively. It is worth noting that a square autocorrelation window with dimensions of 10.7 mm and a window step of 1.5 mm were used in these measurements.

The results demonstrate that the simple autocorrelation in the x-direction fails to differentiate between lesions within the brain phantom, both for the X-motion displacement field (see Figure 11a) and the Z-motion displacement field (see Figure 11d). The application of simple autocorrelation in the y-direction on the X-motion dataset (see Figure 11b) leads to the visualization of two lesions in the background material. However, in the Z-motion dataset (see Figure 11e), only one of the lesions is visualized, along with two artificial unrealistic lesions. In contrast, as is evident in Figure 11c,f, the estimated SWS maps using AIA for the X-motion and Z-motion displacement fields successfully visualize the presence of two lesions within the brain-mimicking phantom. Furthermore, the left-side lesion (18 mm diameter) appears as a larger zone compared to the right-side lesion (12 mm diameter), indicating that the left-side lesion is larger in size than the right-side lesion.

Table 7. SWS in the background material for the brain phantom MRE estimated using simple autocorrelation approaches and the AIA approach in X-motion and Z-motion displacement fields, errors were calculated using the ground truth SWS value of 1.83 (m/s) for the brain phantom as reported in Kabir et al. [47].

Displacement Field	SWS (m/s) in Background Using Simple Autocorrelation (kx)	SWS (m/s) in Background Using Simple Autocorrelation (ky)	SWS (m/s) in Background Using AIA	SWS Estimation Error Using Simple Autocorrelation (kx)	SWS Estimation Error Using Simple Autocorrelation (ky)	SWS Estimation Error Using AIA
X-motion	4.80 ± 1.22	1.19 ± 0.13	1.76 ± 0.15	162%	35%	4%
Z-motion	4.59 ± 1.44	2.09 ± 0.48	1.83 ± 0.22	151%	14%	0%

presents the SWS for the MRE brain phantom in the background material, estimated using simple autocorrelation approaches and the AIA approach in X-motion and Z-motion displacement fields. All MRE measurements of SWS and SNR were conducted considering the almost-full 2D background, with some distance to the lesions’ boundaries. According to Kabir et al. [47], the background shear modulus for the brain phantom is 3.34 ± 0.04 kPa, equivalent to SWSs of 1.83 m/s. It is evident that the estimated SWS using kx is an outlier with an average error of 157% for the X-motion and Z-motion displacement fields in the background. This error is 25% using ky. In contrast, the average error for estimating SWS using AIA in the background is 2%.

Table 8. SNR of estimated SWS in the background material for the brain phantom MRE using simple autocorrelation approaches and the AIA approach in X-motion and Z-motion displacement fields.

Displacement Field	SNR (dB) in Background Using Simple Autocorrelation (kx)	SNR (dB) in Background Using Simple Autocorrelation (ky)	SNR (dB) in Background Using AIA	AIA Improvement over kx	AIA Improvement over ky
X-motion	5.96	9.71	10.66	79%	10%
Z-motion	5.03	6.43	9.15	82%	42%

illustrates the SNR of estimated SWS within the background material for the brain phantom MRE using simple autocorrelation approaches in the x and y directions, as well as the AIA approach for both X-motion and Z-motion displacement fields. These results clearly demonstrate a substantial enhancement in using the AIA approach compared to simple autocorrelation approaches.

4.4. Optical Coherence Elastography

In the OCE experiment, a 3D dataset was obtained depicting the displacement field in the optical direction of the system, denoted as the z-direction in this study. The SWS maps for the phantom including a stiff lesion were determined by employing simple autocorrelation estimators in the x and y directions, along with the AIA approach. In this case study, the autocorrelation window size was adjusted to 1.3 mm × 1.3 mm.

Figure 12 presents the SWS obtained from the simple x-direction autocorrelation (see Figure 12a), simple y-direction autocorrelation (see Figure 12b), and the AIA approach (see Figure 12c) for the phantom featuring a stiff lesion at the frequency of 1500 Hz. The lesion’s location and size are outlined by black dashed lines in the center of the phantom. As shown in Figure 12, the AIA approach effectively visualizes the stiff lesion as a distinct red region within the soft background. In contrast, the simple autocorrelation approaches encountered challenges in accurately identifying the stiff lesion. Additionally, simple autocorrelation approaches resulted in a significant level of uncertainty in SWS estimations within the background region.

Table 9. SWS in the background material the gelatin phantom OCE estimated using simple autocorrelation approaches and the AIA approach, errors were calculated using the ground truth SWS value of 1.94 (m/s) for a 5% gelatin phantom as reported in Zvietcovich et al. (2019) [35].

Excitation Frequency	SWS (m/s) in Background Using Simple Autocorrelation (kx)	SWS (m/s) in Background Using Simple Autocorrelation (ky)	SWS (m/s) in Background Using AIA	SWS Estimation Error Using Simple Autocorrelation (kx)	SWS Estimation Error Using Simple Autocorrelation (ky)	SWS Estimation Error Using AIA
1500 Hz	2.43	2.28	2.05	25%	18%	6%

displays the comparison of SWS within the background material of the gelatin phantom estimated using simple autocorrelation techniques and the AIA approach. The estimated SWS and SNR in OCE analysis were computed considering the almost-full 3D background with some distance to the boundaries. Errors were computed using the established ground truth value of 1.94 (m/s) for the SWS of the 5% gelatin phantom, as documented in Zvietcovich et al. [35].

Table 10. SNR of estimated SWS in the background material for the gelatin phantom OCE using simple autocorrelation approaches and the AIA approach.

Excitation Frequency	SNR (dB) in Background Using Simple Autocorrelation (kx)	SNR (dB) in Background Using Simple Autocorrelation (ky)	SNR (dB) in Background Using AIA	AIA Improvement over kx	AIA Improvement over ky
1500 Hz	4.91	4.69	10.66	117%	127%

displays the SNR of estimated SWS in the background material for the gelatin phantom OCE at the frequency of 1500 Hz utilizing simple autocorrelation approaches and the AIA approach. This table highlights the substantial improvement offered by AIA over simple autocorrelation approaches, with a 117% improvement over kx and a 127% improvement over ky. These results highlight the exceptional capabilities of the AIA approach in effectively identifying and characterizing stiff lesions within the gelatin phantom in OCE experiments.

5. Conclusions

This study introduced the AIA estimator, developed explicitly for elastography measurements in reverberant shear-wave fields. Through the integration of data from all directions of the 2D autocorrelation function, the AIA method exhibited robustness in the presence of data noise. We note that the 1D AIA approach is more computationally efficient than the alternative 2D curve fitting to the theoretical autocorrelation function in any imaging plane, where 2D computations increase as N² points within the window (as opposed to N points for AIA). Furthermore, the N² fit has a greater potential for local minima within the mean-squared error profile. Finally, as shown in Figure 8b, the autocorrelation function in imperfect reverberant fields (having a dominant direction) can have an anomalous pattern in one direction that does not readily converge to the theoretical curve, and this is ameliorated in the AIA approach.

A comprehensive examination was conducted to assess the effectiveness of this advanced autocorrelation approach, using a numerical simulation featuring a stiff branching tube in a uniform background. The AIA approach demonstrated enhanced accuracy in estimating the SWS ratio between the stiff branching tube and the background material. Furthermore, the SNR of the estimated SWS within the background material illustrates the improvement achieved with the AIA approach compared to simple autocorrelation approaches. The practical performance of the AIA approach was assessed through a series of experiments, including ultrasound elastography of a breast phantom with a lesion, ultrasound elastography of the liver–kidney region, MRE of a brain-mimicking phantom with two lesions, and OCE of a gelatin phantom with a stiff lesion. Across this diverse array of experiments, shear-wave fields, and excitation frequencies, the AIA results consistently demonstrated substantial enhancements in the estimated SWS as well as the SNR values. This study has presented the AIA estimator as a robust tool for improving SWS estimation in SWE, offering enhanced accuracy across a spectrum of applications, even in the presence of imperfect reverberant shear-wave fields. These findings highlighted the potential of the AIA approach to advance the field of elastography and contribute to a more precise characterization of the elastic properties of tissue in clinical applications. However, due to the integration of data in all directions, this approach will not resolve the axes of anisotropic materials with differing material properties in different directions, and, in those cases, we would switch the analysis of reverberant shear-wave fields to the approach of Aleman et al. [40].