Phase-Aberration Correction in Shear-Wave Elastography Imaging Using Local Speed-of-Sound Adaptive Beamforming

Shear wave elasticity imaging (SWEI) is a non-invasive imaging modality that provides tissue elasticity information by measuring the travelling speed of an induced shear-wave. It is commercially available on clinical ultrasound scanners and popularly used in the diagnosis and staging of liver disease and breast cancer. In conventional SWEI methods, a sequence of acoustic radiation force (ARF) pushes are used for inducing a shear-wave, which is tracked using high frame-rate multi-angle plane wave imaging (MA-PWI) to estimate the shear-wave speed (SWS). Conventionally, these plane waves are beamformed using a constant speed-of-sound (SoS), assuming an a-priori known and homogeneous tissue medium. However, soft tissues are inhomogeneous, with intrinsic SoS variations. In this work, we study the SoS effects and inhomogeneities on SWS estimation, using simulation and phantoms experiments with porcine muscle as an abbarator, and show how these aberrations can be corrected using local speed-of-sound adaptive beamforming. For shear-wave tracking, we compare standard beamform with spatially constant SoS values to software beamforming with locally varying SoS maps. We show that, given SoS aberrations, traditional beamforming using a constant SoS, regardless of the utilized SoS value, introduces a substantial bias in the resulting SWS estimations. Average SWS estimation disparity for the same material was observed over 4.3 times worse when a constant SoS value is used compared to that when a known SoS map is used for beamforming. Such biases are shown to be corrected by using a local SoS map in beamforming, indicating the importance of and the need for local SoS reconstruction techniques.

1 Introduction

Shear wave elasticity imaging (SWEI) is a non-invasive imaging technique that maps shear-wave speed (SWS) in tissues. Conventionally SWEI is performed in two steps: First in the vicinity of soft tissue to be imaged, a remote “push” is generated using acoustic radiation force (ARF) to induce shear-waves. Second, these shear-waves are observed using ultrasound imaging to capture lateral shear-wave travel speed [1], to relate this to the underlying tissue shear modulus. Supersonic shear-wave imaging (SSI) aims to estimate shear-wave map of the soft tissue with high signal-to-noise ratio (SNR), by generating shear-waves with the constructive interference of multiple pushes along the depth direction, while tracking the shear-waves using ultrafast plane wave imaging (PWI) [2] typically at 10,000 frames per second. SWS measurements have been used in many clinical applications including the diagnosis and staging of diseases in the liver, breast, and kidney [3]; [4]; [5]. SWEI methods usually assume soft tissues are acoustically homogeneous with a nearly constant SoS, for both the generation of ARF pushes and beamforming of PWI in spatially tracking shear-waves. However, soft tissues are acoustically inhomogeneous, which may thus introduce artifacts in SWS estimation. Shi et al. [6] compared SWS and SWS dispersion values measured in-vivo and ex-vivo on three porcine livers to study the confounding effects of porcine skin/fat/muscle on SWS measurements. Carrascal et al. [7] studied phase aberration and ultrasound attenuation effects on SWS and shear-wave frequency domain characteristics. Huang at al. [8] reported that SWS estimation errors due to phase aberrations originate mainly from tracking rather than the ARF push generation. In this work, we aim to correct such errors in SWS estimation. Several data acquisition sequences have been proposed in the literature to mitigate phase aberration effects in shear-wave tracking [9]; [10]; [11]. However, these methods target mediums with slight variations in SoS, whereas in a clinical setting, several layers with largely varying thicknesses and SoS may exist between the ultrasound transducer and the location of measurement. To alleviate large phase aberration effects, knowing the local SoS distribution would be essential.

Several methods have been introduced in the literature to estimate local SoS distributions in soft tissues, mostly aiming for diagnostic purposes that may be afforded by SoS contrast. These methods are known as ultrasound computed tomography (USCT) methods. Conventional USCT systems are based on submerging the target anatomical structure in a water bath, which is equipped with a large number of cylindrically/spherically positioned transducer elements at known locations [12]; [13]. Such transmission USCT systems have great potential for in-vivo breast cancer screening. So it is naturally beneficial to develop SoS imaging to be compatible with existing conventional ultrasound transducers in order to avail several logistic advantages of commercial transducer arrays, also for SoS imaging in the clinics. Time-of-flight recordings together with a passive acoustic reflector [14]; [15] or minute misalignments between images viewed from different angles [16]; [17]; [18] were used for tomographic reconstruction of SoS. Given SoS maps, delays to any spatial location can also be calculated to correct for aberrations caused by SoS inhomogeneities; these delays can be used for beamforming, called SoS-adaptive beamforming, which was shown to increase the resolution of B-mode imaging [19]. Herein, we hypothesize that SoS-adaptive beamforming may be used for mitigating effects on SWS estimation, e.g., via hindering displacement estimation in shear-wave tracking and/or confounding the apparent speed of shear wavefront in consecutive frames due to aberrations. We demonstrate this experimentally using simulations and ex-vivo phantoms that are designed to introduce aberration effects.

2 Methods

2.1 Data Acquisition

In this study, we employ shear-waves induced via ARF using the supersonic shear-wave imaging technique [2], which generates a quasi-cylindrical shear-wave front. We use five consecutive high-intensity ARF pushes of each 200  μs duration, at five axially separated foci with a separation of 2.5 mm in depth, as illustrated in Figure 1A. Laterally propagating shear-wave is then tracked utilizing multi-angle plane-wave imaging (MA-PWI) at a high frame-rate of 10 K frames per second, using three angled plane waves at (−8°, 0°, 8°) as depicted in Figure 1B. Both the acquisition sequences for generating ARF and for tracking shear-wave imaging using MA-PWI were programmed in a research ultrasound machine (Verasonics, Seattle, WA, United States) with a 128-element linear-array transducer (Philips, ATL L7-4) operated at 5 MHz center frequency.

FIGURE 1

FIGURE 1. Schematic of shear-wave elasticity imaging: (A) shear-wave generation; (B) tracking the shear-wave using multi-angle plane wave imaging (MA-PWI); and (C) data processing to estimate the shear-wave speed (SWS).

2.2 Processing Pipeline

We follow the SWEI data processing pipeline as illustrated in Figure 1C. The raw radio frequency (RF) signals for the MA-PWI frames were collected and recorded during the experimental runtime, in order to software-beamform these using alternative approaches retrospectively. After beamforming, a moving window of three PWs (with distinct angles) is coherently compounded to increase SNR [9]. Note that this does not reduce the frame-rate. From the compounded frames, axial displacements are estimated using a 2D Loupas autocorrelation method [20]. These are then directionally filtered to separate the left and right propagating waves [21] along the depth axis z. Using these shear-wave axial displacement profiles, there are several methods in the literature to estimate the SWS [22]; [23]; [24]; [25]. In this study, we used the 1D SWS estimation method of [25], which we observed to be a robust estimator during preliminary tests for our earlier work [26]. In this method, SWS at a position (x, z) in axial and lateral axes, respectively, is computed using a model of wave propagation, i.e., via multiple normalized cross-correlations (NCC) between w consecutive displacement profile pairs, each p-pixels apart, within a lateral window $\left[x+\frac{w}{2},x-\frac{w}{2}\right]$ at depth z. The final SWS is then the NCC-weighted average of maximum correlating profile delays, normalized by the wave travel distance equivalent to p-pixels. For a detailed implementation of this method, please refer to [25].

Parameter w provides an averaging effect, increasing SNR but hampering resolution by smoothing out spatial variations as a tradeoff. Setting parameter p is also a tradeoff; small values allow comparing displacement waveforms of sufficient similarity (i.e., reducing dissimilarity due to dispersion, etc.), while large values allow increasing precision thanks to waveforms with larger shifts in time. Parameters p is set based on the minimum expected SWS in the given medium, which we assumed to be 1.5 m/s in our study. In this study, we accordingly set w = 4 and p = 8.

2.3 Beamforming

For the beamforming of MA-PWI seen above, we employ delay-and-sum using dynamic aperture with an f-number of 0.75 and using the time-delays computed by one of the following two options: First, as the conventional approach, a constant SoS (e.g., an average value from the literature for the imaged anatomy) is used for computing time delays. When the target anatomy is unknown, a typical value of choice is 1,540 m/s, which is also the SoS for the CIRS phantom used. Given the aberrator and different compositions in imaging field, we tested several constant background values in this work. Second, as an alternative approach, we used known local SoS maps, annotated manually from the images, for aberration correction in beamforming as in [19]. Delays τ_i from all 128 transducer elements to all the locations on an N_x × N_z beamforming grid can be computed given local SoS σ_map using the relationship

${\tau }_i=L_i{\sigma }_{map},i=[\mathrm{1,2},\dots ,N_x\times N_z\times 128],(1)$

where L_i are the rows of a path matrix as in Rau et al. [19]. Note that such path matrix for beamforming only needs to be computed once given a transducer geometry and beamforming grid (depth), so it can also be precomputed. Furthermore, given the fixed image composition across all MA-PWI frames, time delays τ_i are fixed among these\ frames.

3 Experiments

For experiments, we used a standard elasticity phantom, CIRS Elasticity QA (Norfolk, VA, United States), with a manufacturer-declared SoS of 1,540 m/s. To mimic aberrations from real tissues, we placed an ex-vivo porcine muscle on the CIRS phantom and placed the US transducer above it so that the aberration source is in one half of the US imaging field-of-view, as seen in Figure 2A. We filled water on top of the phantom as acoustic coupling medium. Before starting experiments, we ensured that the temperature of the water and hence the muscle sample within are stabilized at 22.4°C, tracked using a thermometer for control.

FIGURE 2

FIGURE 2. Illustration of experimental configurations. B-mode images of water and porcine muscle layers placed on CIRS phantom are shown in (A) for the experimental configurations ① and ②, and in (B) for ③ and ④. Acoustic radiation force push markings are overlaid. Manually-segmented custom SoS maps (σ_map) are shown in (C) for ① and ②, and in (D) for ③ and ④.

In total, we used four experimental settings, with two phantom configurations, as shown in Figure 2. Utilizing experiments with the same material but with different combinations allowed us to study beamforming with different SoS effects on SWS estimation, while differentiating effects on the imaging stage from the ARF push path. We created two phantom configurations by changing the location, volume, and cross-sectional profile of the porcine muscle sample overlaid on the CIRS phantom, as shown in Figures 2A,B. The two configurations show some variation in acoustic propagation characteristics and thereby help us identify any systematic effects and demonstrate repeatability of any findings. For each configuration, we then conducted SWS estimation experiments with ARF pushes (i) using the first 64 elements of the transducer over the water layer side and (ii) using the last 64 elements over the muscle sample side. In both cases, all 128 elements were used for imaging, i.e., to receive the signals during MA-PWI (although not all elements contribute to all beamformed points due to dynamic focusing).

To be used in aberration correction with local SoS maps, we manually segmented the three material regions from the images of the two phantom settings, as seen in Figures 2C,D. We set their SoS as follows: For the CIRS phantom, we used the declared value. For the water, we computed the SoS given our monitored temperature using [27]. For the muscle sample, we measured the time-of-flight (ToF) from the strong reflection at the phantom-muscle interface (see the white lines within markings and ④ in Figures 2A,B) and observed its shift with and without the muscle sample being there, to get a relative estimate of the muscle (given its observed thickness) with respect to the known water SoS. Accordingly, the water, CIRS phantom, and the muscle sample were set to have SoS values of 1,490, 1,540, and 1,580 m/s, respectively. For each acquisition, to generate a plane-wave in 2D (quasi-conical wave in 3D), a total of 5 consecutive ARF pushes were conducted at depths {20.0, 22.5, 25.0, 27.5, and 30.0} mm. In our preliminary experiments with homogeneous medium, SWS estimation was observed to be minimally affected by the SoS choice in ARF push Tx delays. Therefore, we used a fixed SoS of 1,540 m/s for all ARF push Tx delay calculations. SWS was estimated using [25] after beamforming the MA-PWI using constant SoS values of {1,480, 1,490, …, 1,590} m/s as well as using the local SoS map [19]. For each of the four experimental configurations and for any given beamforming setting, the acquisition process described above was repeated five times and the resulting five SWS maps were pointwise averaged to increase the SNR. These averaged SWS maps are used below to compare and quantify effects from SoS aberrations.

4 Results and Discussion

In the following, we first show and analyze the results from the experimental setting in detail, and then summarize the findings from all four settings comparatively.

4.1 Detailed Analysis of a Sample Experimental Setting

To illustrate the effect of aberration when beamformed with different SoS values, the axial velocity profiles at a selected depth of 25 mm (as the middle push location, with a relatively ideal shear wavefront), averaged over an axial window of 3 mm around the selected depth 25 mm, are shown for the experimental configuration in Figure 3(I). These profiles are shown as a function of time across the phantom width, and they were computed from the estimated displacements, for the beamforming of which we first used different constant SoS values (Figures 3(I,a–f)). It is seen that during the propagation of the shear-wave, it becomes distorted along a vertical band around the mid-line of the phantom (x = 0), irrespective of the constant beamforming SoS value. Such distortion band is where part of the received echoes, used in the beamforming, goes through the muscle, while the other half goes through the water, with substantially different SoS values; therefore, the beamforming is largely out-of-phase. Indeed, the severity of the distortion does not seem to change given the beamforming SoS value and even for the manufacturer-declared 1,540 m/s; such distortion is observed as shown circled in Figures 3(I,d). To further illustrate this effect, ribbon plots in Figure 3(II) zoom in on such distortion regions (delineated with a rectangle in Figures 3(I,b)) to visualize the time-domain characteristics of the axial velocity profiles; that is, the observed shape and time shifts between consecutive axial velocity profiles. One can see in these profiles that the velocity profiles in the distortion bands lose the coherent wavefront appearance that they otherwise had.

FIGURE 3

FIGURE 3. Experimental results. (i) Shear-wave axial velocity time and propagation distance profiles; (ii) Axial velocity profiles corresponding to rectangular region (i,b); (iii) resulting shear-wave speed (SWS) maps; and (iv) SWS profiles over selected rectangular regions of interest (ROI), with the errorbars representing the standard deviations of SWS values along the vertical axis, with the dashed lines marking the average of their means, black for the left side and green for the right. Each column shows results for beamforming the same data acquisition with different speed-of-sound (SoS) maps (A–F) for constant SoS maps from 1,480 to 1,580 m/s, and (G) using a spatially varying SoS map σ_map. We do not report the results, by masking them in white, within vertical bands around the lateral push locations where the shear-wave travels before the tracking starts.

Figure 3(III) displays the resulting SWS maps obtained within a rectangular field-of-view (FoV) of size 14 × 35.2 mm ([22,36] mm axially and [−17.6,17.6] mm laterally). Such an FoV was selected where we know that the (CIRS) phantom has a constant homogeneous SWS, given the phantom construct and the survey B-mode image in Figure 2. In the estimated SWS maps, severe distortions can be observed around the mid vertical lines (x = 0) due to the large SoS inhomogeneity; i.e., large variations in backscatter wavefront paths are considered in beamforming, which introduce aberrations in the shear-wave particle velocity profiles, regardless of the depth where the shear wavefront is observed. Therefore, SWS cannot be reliably estimated along such narrow vertical bands. Besides this immeasurable region, a striking observation is that the SWS estimated on the left and that on the right sides of such band are not equal (seen as different blue hues), although the estimations can be made with relatively high SNR. This is also observed in the velocity profiles in Figure 3(II) where the wavefronts on either side of the distortion band exhibit a mismatch in gradient (which leads the local SWS estimate), as illustrated by marking their automatically selected peak locations in Figures 3(II,d). This is not only a random noise or distortion and indicates a systematic bias in SWS estimates due to an aberrator and the conventional beamforming process, which may affect and alter any diagnostic decision in a clinical setting. Throughout the rest of this article, we systematically analyze and further quantify such effects and study if such bias can be reduced or removed.

To better illustrate the differing SWS values due to aberration, Figure 3(IV) presents some statistics of SWS values within rectangular ROIs of 3.0 × 4.5 mm that are placed equidistant from the middle of the imaging width, as marked in Figures 3(III,b). In these subfigures, the values around the distortion bands were not plotted as these explode the axis range and also have large very variations due to low estimation accuracy. From the SWS mean values, marked with dashed lines of black and green in Figure 3(IV), one can see that the SWS estimations on either side differ largely, despite both sides having relatively low standard deviations. Since we know from the phantom construct that the given FoV has a constant homogeneous SWS, one would expect for the SWS values to match. Note that even if the beamforming SoS errors would change the arrival time of echos to the transducer (and therefore potentially degrade image quality and resolution), one would not expect this to affect the speed estimation of laterally propagating shear-waves and certainly not in a systematic or an easily explainable way. Further of note is that this systematic bias occurs regardless of the chosen constant SoS value for beamforming, seen across the columns (a–f). To the best of our knowledge, we are unaware of a systematic reporting or study of this in the literature.

For beamforming the MA-PWI data using a spatially varying SoS map in Figure 3(g), we observe in (I,g) and IIIg that the distortion regions are substantially minimized. From the continuous appearance of the velocity profiles (peaks) in (II,g) as well as the constant color hue in the SWS map in (III,g), we can see that using the knowledge of the local SoS variation in beamforming can indeed remove the SWS estimation bias. This is further illustrated in (IV,g) where SWS means within ROIs on both sides have statistically indifferent SWS estimations, showing that beamforming with local SoS maps can remedy such SWS bias and problem.

4.2 Comparison of Beamforming for All Experimental Settings

To validate the above observations by controlling for the medium that the ARF push beams have travelled as well as any systematic effect from the direction of the phantom/muscle placement, we analyzed the results from four different experimental configurations. We tested two phantom configurations, by changing the muscle to the left (in the experiments ① and ②) and the right side (in ③ and ④). In each configuration, we applied the ARF push either through the water layer (in ① and ③) or through the porcine muscle (in ② and ④). Since the configurations were nearly symmetric, we used the same ROI boxes above, fixed for all the experimental configurations. We then report the mean and standard deviations of SWS values within each ROI and over five repetitions in Table 1. Note that given the experimental configuration, an ROI is either nearer or farther from the ARF push location, and an ROI is imaged either mostly through water or through muscle. To help isolate any effects from different factors, we use the following terminology in reporting the results below. We denote the SWS results in ROI nearer and farther from the ARF push with S and S′. For results in ROI imaged through water (W) and muscle (M), we use the relevant letter in the subscript, i.e., S_W and S_M. For instance, in Figures 3(IV,d), $S_{\text{W}}>S_{\text{M}}^{\prime }$ although in (IV,d) with local SoS map $S_{\text{W}}\approx S_{\text{M}}^{\prime }$, we tabulate the SWS estimations for different beamforming SoS values in Table 1.

TABLE 1

TABLE 1. Measured SWS values in ROI nearer (S) and farther (S′) from the ARF when imaged through water (i.e., S_W) and muscle (i.e., S_M). Results are shown after beamforming with difference constant SoS as well as a locally-varying SoS map (σ_map). Minimum disparity D achieved across the beamforming settings is highlighted in bold for each experiment.

From the table, it is seen that the SWS estimations are lower in ROIs on the muscle side of the distortion band, in comparison to the water side. This is true irrespective of the ARF push going through the muscle or the water and again irrespective of the muscle side being closer to or farther from the ARF push. We also report the mean absolute difference (disparity) D = |S − S′| between the two ROIs in each experimental setting and SoS. It is seen that for constant SoS values for beamforming, SWS differences D are quite large. For visual comparison of these values, we have plotted SWS values obtained from the mentioned ROIs for each experimental configuration in Figure 4 . Ideally, this SWS disparity should be near zero, since it is the exact same material being measured in both ROIs. Indeed, the average disparity of all four experiments for a typical (also CIRS-declared) SoS of 1,540 m/s is 0.52  m/s, which is over 4.3 times higher than the average disparity (0.12  m/s) when a local SoS map σ_map is used for beamforming the MA-PWI data. Furthermore, even for using a best constant SoS value to minimize such difference for each individual experiment (underlined in the table), the average is 0.42  m/s, which is again 3.5 times higher than that by using the local SoS map. Even if one assumes that a global SoS estimator may be used to predict such constant SoS value, it is seen that this “best” SoS value is not consistent across different experiments; thus, no constant value would be ideal.

FIGURE 4

FIGURE 4. Visual summary of Table 1 with SWS values obtained from the two ROIs in each experimental setting when beamformed with constant SoS values and σ_map.

It is seen that in all the experimental configurations, the mean absolute difference D is the least when MA-PWI data is beamformed using a local SoS map σ_map, with the reported differences being below the reported standard deviations, i.e., the experimental noise floor. It is also worth to note that the SWS values estimated using σ_map are relatively consistent across the experimental setups as well, although these vary largely for the values estimated using constant SoS beamforming. In addition, from the reported standard deviations, it is observed that SWS values exhibit more variations closer to the ARF push, with similar observations reported earlier by us in [28]; [26] as well as by other groups in [29]; [30].

In this work, we corrected the phase aberrations in the received MA-PWI data but not in transit delays for correcting phase aberrations of ARF pushes. Therefore the observed biases in SWS estimation are due to the cumulative effect of both phase aberrations in generation and tracking shear-waves. Nevertheless, in Table 1, the SWS values reported in the CIRS phantom beneath the water and muscle sides are seen to be consistent between experiments ① and ②, despite the ARF pushes having used fixed Tx delays while being applied through mediums with different SoS. Similar observations can be made from experiments and ④, indicating that phase aberrations in ARF push generation potentially affect the generated and observed SWS minimally. These observations also corroborate the findings in [8] as well as our preliminary experiments where incorrect SoS values in ARF generation did not result in significant differences in SWS estimation. It is however unclear if such aberrations may cause differences in other SWE applications, such as frequency-dependent, non-linearity, or attenuation measurements, e.g., by changing the spectral content of induced shear-waves. Nevertheless, one could also correct ARF Tx delays using local SoS maps σ_map, if this becomes relevant for a given application scenario. Aberration-based SWS errors observed for tracking with PWI are expected to also be present if focused beams are used, since PW and focused beams travel through similar tissue regions, especially in a layered medium. These SWS errors should similarly be correctable using local SoS adaptive beamforming. Note that in our beamforming of MA-PWI, we are already utilizing Rx-focusing, where the adaptive beamforming is shown to be advantageous, so given reciprocity with time-reversal, one would expect such benefit to also exist for focusing on the Tx side.

To qualitatively illustrate the effect of aberration on shear wave amplitudes, when ARF is focused through different materials and MA-PWI is beamformed with different SoS values, the mean maximum amplitude of axial velocity profiles within the same ROI that was used for plotting axial velocity profiles in Figures 3(I) is shown in Figure 5 for all the experimental configurations. The average of maximum shear-wave particle velocities in experiments ② and ④ with ARF pushes through the muscle are seen to be ≈3.8 times smaller compared to that of experiments ① and ④ with pushes through water, given that ARF focusing was performed with a fixed assumed SoS. These findings corroborate the observations in [7]. In contrast, differences in shear-wave particle velocities due to beamforming with different SoS values are seen to be minimal for all experimental configurations.

FIGURE 5

FIGURE 5. Maximum shear-wave particle axial velocity values for the four experimental configurations, with each line showing the velocity observation from the same MA-PWI data but using different SoS assumptions in beamforming, that is, constant values from 1,480 to 1,580 m/s vs. a spatial map σ_map.

We herein show that local SoS-based adaptive beamforming corrects phase aberration effects on SWS estimation. To that end, robust and accurate local SoS estimation is key. Several groups showed SoS reconstruction for submersible body parts using ring or rotating transducer setups. To alleviate the need for complex setups, hand-held solutions with acoustic reflector based, e.g., [15]; [31], and pulse-echo disparity based, e.g., [17]; [18], tomographic reconstruction methods were also demonstrated. To achieve fast, real-time SoS map estimations, such reconstructions have also been accelerated using deep-learning based techniques in [32] and [33].

4.3 Simulation Study

To better reason about the observations above, we conduct a simulation experiment. For a simulation of deformations caused by shear-waves and the imaging thereof, one would require a complex continuum mechanics simulation that can faithfully model the ARF, tissue dynamics, and the imaging interactions, each of which requires very complex and mostly unknown parameterizations. Not only setting up such a simulation would be very difficult, but also any results would be highly sensitive to chosen parametrizations and thus potentially questionable. Therefore, we chose herein a simplified simulation scheme, where we neglect the ARF mechanics and energy exchange, as well as the shear-wave dynamics, and focus instead on the effects from SoS in beamforming. This is also in line with our observations in Figure 5 that although the maximum particle velocities and hence the displacement amplitudes may vary, mainly due to ARF push mechanics, the estimated SWS nevertheless is relatively independent of such amplitudes and thus ARF mechanics and is largely dependent on the beamforming (SoS) as seen from the rest of our results. Accordingly, we focus on an already-generated shear-wave and simulate the SWS estimation process instead. Most forms of SWS estimation operates on tracked displacement-time profiles (DTPs) as an observation of the passing shear wavefront, as illustrated in Figure 3(II). Methods may use cross-correlation of such DTPs at known spatial increments (variable p in our method), or even simple peak-picking in DTPs (e.g., the black dots in Figures 3(II,D)), in order to infer the time elapsed for the wave to travel this distance. In our simulations, we place scatterers (“anchors”) to mark the physical locations of regular intervals where the shear-wave generates such DTPs. We then look at the beamformed images to see where the anchors (scatterers) are observed, which will also be the location where the DTPs will be recorded from an SWS-estimating observer point-of-view.

We use the SoS map σ_map in Figure 2C in the k-wave acoustics simulation toolbox [34]. To mark physical tissue anchor points along a laterally propagating wavefront with a constant speed, we place a simulated point scatterer at a fixed depth of z = 27.3  mm, while changing its lateral location from x = −13.8 mm to + 13.8 mm, at intervals of 0.3 mm, leading to 93 individual scatterer locations. For each scatterer location, an MA-PWI simulated with three angles {−8, 0, 8} was acquired, with each angle beamformed using the delay-and-sum method, assuming a homogeneous medium with constant SoS of 1,580 m/s, and these triplet images were then spatially compounded as in our standard SWS imaging approach above. In each of the 93 beamformed and compounded frames, we then easily identify the peak of the isolated scatterer (based on envelope intensity and a quadratic subsample approximation [35]) and treat this apparent anchor location as the observed location of a DTP which actually occurs in the physical location of the anchor, as seen in Figure 6. Assuming an MA-PWI framerate of 10 K frames/s, 0.3 mm physical anchor distance per tracking frame time of 100 μs corresponds to a 3 m/s shear-wave travel.

FIGURE 6

FIGURE 6. Illustration of simulation. A given SWS encountered as displacement-time profiles (DTPs) at physical anchor (scatterer) location (left) may appear as a different SWS if the same DTPs are observed at different image locations (right), e.g., due to aberrations in beamforming.

In Figure 7A, we plot the true axial position given in the simulation against the apparent axial-axis position measured from the beamformed data. As can be seen, underneath the muscle tissue (the SoS of which is closer to the chosen constant beamforming SoS), the apparent depth of the scatterer is closer to the actual depth. Note that a near-vertical planar shear wavefront, such axial shift would not have a major effect. To gain insight into lateral propagation characteristics, we also plot the true lateral positions against the apparent lateral apparent position in Figure 7B. One can see there the offset from the true lateral positions in a large transition region where some but not all beamforming paths cross the aberrator. This then perturbs the SWS observed from these observations. An extreme example, for instance, the lateral observations of two consecutive locations x = {−2.25, −1.95} mm that are then observed as the apparent location x = {−1.35, 0.75} mm, which using some finite-difference speed estimation, hence would distort the actual shear-wave speed of 3 m/s, to locally appear as 21 m/s in the beamformed data. For SWS estimation in practice, instead of finite-differences, often some finite filter length is used for noise suppression. In Figure 7C, we converted the lateral position observations to SWS estimation using a smoothing approach similar to the practical phantom experiment earlier. This figure illustrates how SWS is inflated or deflated due to SoS inhomogeneity effects although the actual physical SWS is 3 m/s. In this figure, we mark the regions corresponding to the width of the left and right ROIs reported earlier in our experiments, with the marked lines indicating the mean SWS estimation in these regions. From this simulation result, a slight overestimation of SWS under the water side as well as a gross underestimation under the muscle side is observed. These results fully corroborate the findings in the experiments in Figure 3(IV).

FIGURE 7

FIGURE 7. Simulation results. True vs. apparent axial (A) and lateral (B) positions of a scatterer, observed via beamforming with a homogeneous SoS map of 1,580 m/s. Apparent lateral speed (C) computed using the estimators similarly to the experimental case. The black markings represent the mean estimated SWS within the width of the two ROIs used in the experiments and delineated in Figure 3(III,b).

Note that if we see the anchors at different locations in the beamformed image, we would also observe a corresponding DTP at those locations in a SWS tracking scenario. Therefore, a cross-correlation yielding the same time delay would be observed at a different apparent distance, hence resulting in differently observed SWS. For instance, if two neighbouring anchors are to be observed laterally 0.5 mm apart, since the cross-correlation of their DTPs would yield a time lapse of 100 μs, we would observe an SWS of 5 m/s, which describes the value observed around x = −2.5 mm in the simulation results in Figure 7. By representing the observed DTP locations with a scatterer, we omit the effects of SoS in the displacement estimation itself. This is based on the assumption that the SoS-related beamforming delays can differ only negligibly across the spatial range of shear-wave related particle displacements. Furthermore, any delay differences that may skew the amplitude of observed particle velocities are likely to affect the SWS estimation negligibly, as also observed from the particle amplitudes reported in Figure 5 being unrelated to our resulting SWS estimations. In our simulations, we ignored DTP distortions due to aberrations, i.e., cross correlation leading to any artifactual time shifts in addition to spatial shifts. Such distortions occur, e.g., beneath refracting edges as seen within the ±1 mm band in the DTPs in Figure 3(II). Nevertheless, such distortions often lead to random noise with non-systematic offset and to very low correlation coefficients with SWS estimation being inconclusive in these distorted regions. These assumptions need to be further studied in the future with a comprehensive full-pipeline simulation study.

5 Conclusion

In this work, SoS aberration effects on SWEI when beamformed with different SoS have been studied, with ARF pushes and shear-wave observations through water and muscle samples. We found that ARF push generation affects the generated as well as observed SWS minimally, while the beamforming errors and aberrations may significantly alter the observed SWS values. Average SWS disparity when a constant SoS value is used for beamforming was found over 4.3 times worse than using a known SoS map. Indeed there was no single constant SoS value that could mitigate the observed SWS biases. Nevertheless, it is shown that using a known SoS map in beamforming delay calculations prevents such SWS estimation biases. Note that even using only direct linear paths in beamforming, that is, not taking potential refractions into account, major SWS errors could be prevented.

Without aberration correction, the observed differences in measured SWS were substantial. In terms of elastic moduli, e.g., computed given the SWS values in Table 1 at 1,540 m/s and assuming density of 1 g/cm3, our results show that for the same material, one may find 7.56 vs. 4.67 kPa for experiment–ARF through W and 8.94 vs. 5.15 kPa for experiment–ARF through M. These are important differences which may cause misdiagnosis, especially given that absolute values of SWS measurements are typically used in staging diseases in the liver, breast, and kidney [3]; [4]; [5]. These large differences indicate the difficulty in standardizing SWS measurements and that diseases may be staged incorrectly which can impact management and treatment decisions. For instance, with 7 kPa being the common threshold between mild (F1) and significant (F2) fibrosis, a patient may be staged differently given the pairs of measurements above. It needs to be further studied if such measurement errors would scale linearly or stay constant at increasing SWS values.

Since soft tissues are intrinsically inhomogeneous, for accurate tissue characterization and diagnosis using SWEI, it appears imperative given our study to beamform MA-PWI data using accurate SoS distributions of the medium to alleviate possible confounding effects of SoS and beamforming on estimated SWS. For estimating such local SoS distributions, one could use the same ultrasound transducer used for SWEI as in [16]; [17]; [18] as demonstrated for aberration correction in [19].

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

BC contributed with conducting experiments, data processing, analysis and discussion of results, and manuscript drafting. RR contributed with supervision, data processing, analysis, and discussion of results. OG contributed with supervision, analysis, discussion of results, and manuscript drafting.

Funding

Funding was provided by the Swiss National Science Foundation (SNSF) grant #179116 and the Promedica Foundation, Chur.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Sarvazyan AP, Rudenko OV, Swanson SD, Fowlkes JB, Emelianov SY. Shear Wave Elasticity Imaging: a New Ultrasonic Technology of Medical Diagnostics. Ultrasound Med Biol (1998) 24:1419–35. doi:10.1016/s0301-5629(98)00110-0

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bercoff J, Tanter M, Fink M. Supersonic Shear Imaging: a New Technique for Soft Tissue Elasticity Mapping. IEEE Trans Ultrason Ferroelect., Freq Contr (2004) 51:396–409. doi:10.1109/tuffc.2004.1295425

CrossRef Full Text | Google Scholar

3. Sarvazyan A, J. Hall T, W. Urban M, Fatemi M, R. Aglyamov S, S. Garra B. An Overview of Elastography-An Emerging Branch of Medical Imaging. Cmir (2011) 7:255–82. doi:10.2174/157340511798038684

CrossRef Full Text | Google Scholar

4. Deffieux T, Gennisson J-L, Bousquet L, Corouge M, Cosconea S, Amroun D, et al. Investigating Liver Stiffness and Viscosity for Fibrosis, Steatosis and Activity Staging Using Shear Wave Elastography. J Hepatol (2015) 62:317–24. doi:10.1016/j.jhep.2014.09.020

CrossRef Full Text | Google Scholar

5. Tanter M, Bercoff J, Athanasiou A, Deffieux T, Gennisson J-L, Montaldo G, et al. Quantitative Assessment of Breast Lesion Viscoelasticity: Initial Clinical Results Using Supersonic Shear Imaging. Ultrasound Med Biol (2008) 34:1373–86. doi:10.1016/j.ultrasmedbio.2008.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Shi Y, Xie H, Shamdasani V, Fraser J, Robert J-L, Zhou S. Phase Aberration in Shear Wave Dispersion Ultrasound Vibrometry. In: 2011 IEEE International Ultrasonics Symposium; Orlando, FL, USA18-21 Oct. 2011. IEEE (2011). p. 2408–11. doi:10.1109/ultsym.2011.0598

CrossRef Full Text | Google Scholar

7. Amador Carrascal C, Aristizabal S, Greenleaf JF, Urban MW. Phase Aberration and Attenuation Effects on Acoustic Radiation Force-Based Shear Wave Generation. IEEE Trans Ultrason Ferroelect., Freq Contr (2016) 63:222–32. doi:10.1109/tuffc.2016.2515366

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Huang S-W, Xie H, Robert J-L, Nguyen M, Shamdasani V. Phase Aberration in Ultrasound Shear Wave Elastography-Impacts on Push and Tracking. In: 2016 IEEE International Ultrasonics Symposium (IUS); 18-21 Sept. 2016; Tours, France. IEEE (2016). p. 1–4. doi:10.1109/ultsym.2016.7728569

CrossRef Full Text | Google Scholar

9. Montaldo G, Tanter M, Bercoff J, Benech N, Fink M. Coherent Plane-Wave Compounding for Very High Frame Rate Ultrasonography and Transient Elastography. IEEE Trans Ultrason Ferroelect., Freq Contr (2009) 56:489–506. doi:10.1109/tuffc.2009.1067

CrossRef Full Text | Google Scholar

10. Amador C, Song P, Meixner DD, Chen S, Urban MW. Improvement of Shear Wave Motion Detection Using Harmonic Imaging in Healthy Human Liver. Ultrasound Med Biol (2016) 42:1031–41. doi:10.1016/j.ultrasmedbio.2015.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Espíndola D, Lee S, Pinton G. Shear Shock Waves Observed in the Brain. Phys Rev Appl (2017) 8:044024. doi:10.1103/physrevapplied.8.044024

CrossRef Full Text | Google Scholar

12. Greenleaf JF, Bahn RC. Clinical Imaging with Transmissive Ultrasonic Computerized Tomography. IEEE Trans Biomed Eng (1981) BME-28:177–85. doi:10.1109/tbme.1981.324789

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Gemmeke H, Ruiter NV. 3d Ultrasound Computer Tomography for Medical Imaging. Nucl Instr Methods Phys Res Section A: Acc Spectrometers, Detectors Associated Equipment (2007) 580:1057–65. doi:10.1016/j.nima.2007.06.116

CrossRef Full Text | Google Scholar

14. Krueger M, Burow V, Hiltawsky K, Ermert H. Limited Angle Ultrasonic Transmission Tomography of the Compressed Female Breast. In: 1998 IEEE Ultrasonics Symposium. Proceedings (Cat. No. 98CH36102); 5-8 Oct. 1998; Sendai, Japan. IEEE (1998). p. 1345–8. doi:10.1109/ULTSYM.1998.765089

CrossRef Full Text | Google Scholar

15. Sanabria SJ, Rominger MB, Goksel O. Speed-of-sound Imaging Based on Reflector Delineation. IEEE Trans Biomed Eng (2018)p. 1949–1962. doi:10.1109/TBME.2018.2881302

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Jaeger M, Held G, Peeters S, Preisser S, Grünig M, Frenz M. Computed Ultrasound Tomography in echo Mode for Imaging Speed of Sound Using Pulse-echo Sonography: Proof of Principle. Ultrasound Med Biol (2015) 41:235–50. doi:10.1016/j.ultrasmedbio.2014.05.019

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Sanabria SJ, Ozkan E, Rominger M, Goksel O. Spatial Domain Reconstruction for Imaging Speed-Of-Sound with Pulse-echo Ultrasound: Simulation and In Vivo Study. Phys Med Biol (2018) 63:215015. doi:10.1088/1361-6560/aae2fb

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Rau R, Schweizer D, Vishnevskiy V, Goksel O. Speed-of-sound Imaging Using Diverging Waves. Int J Compu Assist Radiol (2021):1–11. doi:10.1007/s11548-021-02426-w

CrossRef Full Text | Google Scholar

19. Rau R, Schweizer D, Vishnevskiy V, Goksel O. Ultrasound Aberration Correction Based on Local Speed-Of-Sound Map Estimation. In: 2019 IEEE International Ultrasonics Symposium (IUS). IEEE (2019). p. 2003–6. doi:10.1109/ultsym.2019.8926297

CrossRef Full Text | Google Scholar

20. Loupas T, Powers JT, Gill RW. An Axial Velocity Estimator for Ultrasound Blood Flow Imaging, Based on a Full Evaluation of the Doppler Equation by Means of a Two-Dimensional Autocorrelation Approach. IEEE Trans Ultrason Ferroelect., Freq Contr (1995) 42:672–88. doi:10.1109/58.393110

CrossRef Full Text | Google Scholar

21. Manduca A, Lake DS, Kruse SA, Ehman RL. Spatio-temporal Directional Filtering for Improved Inversion of Mr Elastography Images. Med image Anal (2003) 7:465–73. doi:10.1016/s1361-8415(03)00038-0

PubMed Abstract | CrossRef Full Text | Google Scholar

22. McLaughlin J, Renzi D. Shear Wave Speed Recovery in Transient Elastography and Supersonic Imaging Using Propagating Fronts. Inverse Probl (2006) 22:681–706. doi:10.1088/0266-5611/22/2/018

CrossRef Full Text | Google Scholar

23. Wang MH, Palmeri ML, Rotemberg VM, Rouze NC, Nightingale KR. Improving the Robustness of Time-Of-Flight Based Shear Wave Speed Reconstruction Methods Using Ransac in Human Liver In Vivo. Ultrasound Med Biol (2010) 36:802–13. doi:10.1016/j.ultrasmedbio.2010.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Rouze NC, Wang MH, Palmeri ML, Nightingale KR. Robust Estimation of Time-Of-Flight Shear Wave Speed Using a Radon Sum Transformation. IEEE Trans Ultrason Ferroelect., Freq Contr (2010) 57:2662–70. doi:10.1109/tuffc.2010.1740

CrossRef Full Text | Google Scholar

25. Song P, Manduca A, Zhao H, Urban MW, Greenleaf JF, Chen S. Fast Shear Compounding Using Robust 2-d Shear Wave Speed Calculation and Multi-Directional Filtering. Ultrasound Med Biol (2014) 40:1343–55. doi:10.1016/j.ultrasmedbio.2013.12.026

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Chintada BR, Rau R, Goksel O. Nonlinear Characterization of Tissue Viscoelasticity with Acoustoelastic Attenuation of Shear Waves. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control (2021). doi:10.1109/TUFFC.2021.3105339

CrossRef Full Text | Google Scholar

27. Bilaniuk N, Wong GSK. Speed of Sound in Pure Water as a Function of Temperature. The J Acoust Soc America (1993) 93:1609–12. doi:10.1121/1.406819

CrossRef Full Text | Google Scholar

28. Otesteanu CF, Chintada BR, Rominger M, Sanabria S, Goksel O. Spectral Quantification of Nonlinear Elasticity Using Acousto-Elasticity and Shear-Wave Dispersion. IEEE Trans Ultrason ferroelectrics, frequency Control (2019) 66(12):1845–55. doi:10.1109/tuffc.2019.2933952

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bernard S, Kazemirad S, Cloutier G. A Frequency-Shift Method to Measure Shear-Wave Attenuation in Soft Tissues. IEEE Trans Ultrason ferroelectrics, frequency Control (2016) 64:514–24. doi:10.1109/TUFFC.2016.2634329

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kijanka P, Urban M. Two-point Frequency Shift (2p-Fs) Method for Shear Wave Attenuation Measurement. IEEE Trans Ultrason ferroelectrics, frequency Control (2019) 67:483–96. doi:10.1109/TUFFC.2019.2945620

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Chintada BR, Rau R, Goksel O. Time of Arrival Delineation in echo Traces for Reflection Ultrasound Tomography. In: 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI); 13-16 April 2021; Nice, France. IEEE (2021). p. 1342–5. doi:10.1109/isbi48211.2021.9433846

CrossRef Full Text | Google Scholar

32. Vishnevskiy V, Rau R, Goksel O. Deep Variational Networks with Exponential Weighting for Learning Computed Tomography. In: International Conference on Medical Image Computing and Computer-Assisted Intervention; October 13–17, 2019; Shenzhen, China. Springer (2019). p. 310–8. doi:10.1007/978-3-030-32226-7_35

CrossRef Full Text | Google Scholar

33. Bernhardt M, Vishnevskiy V, Rau R, Goksel O. Training Variational Networks with Multidomain Simulations: Speed-Of-Sound Image Reconstruction. IEEE Trans Ultrason Ferroelect., Freq Contr (2020) 67:2584–94. doi:10.1109/tuffc.2020.3010186

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Treeby BE, Cox BT. K-Wave: Matlab Toolbox for the Simulation and Reconstruction of Photoacoustic Wave fields. J Biomed Opt (2010) 15:021314. doi:10.1117/1.3360308

CrossRef Full Text | Google Scholar

35. Azar RZ, Goksel O, Salcudean SE. Sub-sample Displacement Estimation from Digitized Ultrasound Rf Signals Using Multi-Dimensional Polynomial Fitting of the Cross-Correlation Function. IEEE Trans Ultrason Ferroelect., Freq Contr (2010) 57:2403–20. doi:10.1109/tuffc.2010.1708

CrossRef Full Text | Google Scholar