A PSF-Shape-Based Beamforming Strategy for Robust 2D Motion Estimation in Ultrafast Data

Abstract

This paper presents a framework for motion estimation in ultrafast ultrasound data. It describes a novel approach for determining the sampling grid for ultrafast data based on the system’s point-spread-function (PSF). As a consequence, the cross-correlation functions (CCF) used in the speckle tracking (ST) algorithm will have circular-shaped peaks, which can be interpolated using a 2D interpolation method to estimate subsample displacements. Carotid artery wall motion and parabolic blood flow simulations together with rotating disk experiments using a Verasonics Vantage 256 are used for performance evaluation. Zero-degree plane wave data were acquired using an ATL L5-12 (f_c = 9 MHz) transducer for a range of pulse repetition frequencies (PRFs), resulting in 0–600 µm inter-frame displacements. The proposed methodology was compared to data beamformed on a conventionally spaced grid, combined with the commonly used 1D parabolic interpolation. The PSF-shape-based beamforming grid combined with 2D cubic interpolation showed the most accurate and stable performance with respect to the full range of inter-frame displacements, both for the assessment of blood flow and vessel wall dynamics. The proposed methodology can be used as a protocolled way to beamform ultrafast data and obtain accurate estimates of tissue motion.

1. Introduction

Ultrasound imaging is often used to visualize and estimate motion in our body [1,2,3,4,5,6,7,8,9,10]. Tissue motion and deformation can provide information about the structure, composition and functioning of the tissue and changes in the dynamic behavior of tissue can be used as indicator of disease. For example, stiff regions in breasts are often related to cancers [11,12], the assessment of local arterial deformation provides information about the atherosclerotic progression and vulnerability of lesions [13,14,15,16], increased blood velocities indicate the presence of a stenosis in the vessel [17,18,19], and infarcted myocardial tissue shows changed deformation patterns as compared to healthy tissue [20].

Speckle tracking (ST), first introduced by Trahey et al. [21,22], is one of the motion estimation methods which is frequently utilized. Speckle patterns are tracked to find the displacement of the underlying tissue. Pattern matching techniques are used, where a kernel region in the first acquisition is matched within a surrounding search region in the following acquisition. The location of the best match defines the displacement of the kernel region and thus the displacement of the underlying tissue. ST can be performed using radio frequency (RF) signals [23], envelope signals [24], B-mode images [20,22], or a combination thereof [25]. Bohs et al. [26] wrote a review on the development of this technique, focusing on multi-dimensional flow estimation. They concluded that decorrelation of speckle patterns is one of the major challenges in ST. Strong displacement gradients and velocities with low beam-to-flow angles were described as major sources of speckle decorrelation. Increasing the imaging frame rate can minimize speckle decorrelation [26].

The first published approach for increasing imaging frame rate was by reconstructing multiple image lines simultaneously from each transmit, so-called parallel beamforming [27]. Nowadays, plane wave and diverging wave imaging are used more often, where the transmitted pulse fully illuminates the entire region of interest, enabling full image reconstruction for every transmit event. This allows for imaging at ultrafast frame rates in the kHz range. Hereby, decorrelation of speckle can be minimized and fast moving, sometimes short-lived, motion patterns can be tracked [28,29,30,31,32]. However, frame-to-frame tissue displacements become smaller, since they scale with the imaging frame rate. Consequently, acquiring accurate subsample displacement estimates is of great importance.

The estimation of subsample displacements in ST was already a challenge before the introduction of ultrafast data acquisition, since echo signals are sampled. For conventionally acquired, line-by-line data, interpolation of the pattern matching function is often applied to acquire subsample displacement estimates [25,33,34]. The sampling grid for conventional data is predefined, with normally a line distance equal to the pitch and an axial grid spacing with typically four samples per wavelength. Consequently, the sampling of the pattern matching function is also set. A predefined shape or curve, often parabolic or cosine [34,35], is assumed and fitted to the peak of the function, either in each direction of the 2D function independently [25,35], or by fitting in 2D to obtain a joint estimate of subsample displacement in both directions. Zahiri Azar et al. [33] investigated 2D polynomial fitting to get a joint estimation of subsample displacement and showed significant improvements in terms of bias and standard deviation (SD) compared to the commonly used 1D parabolic and cosine fitting. Interpolating in between echo signals is also often applied to reduce the error of the subsample displacement estimates [25,36,37].

With the use of ultrafast data, subsample displacements estimation remains a challenge. However, it also offers a totally new possibility: since image reconstruction is performed after receiving the raw channel data, the data sampling grid is not fixed and focusing can be performed at each specific location. Redefining the data sampling grid directly influences the sampling of the pattern matching function. Hence, it allows to choose the available samples and their spacing throughout the peak of the pattern matching function, which we will refer to as the appearance or shape of the function throughout this paper. The possibility to influence its appearance allows to optimize the match between the peak of the function and the interpolation method. This will in theory enable more accurate subsample displacement estimation [34].

More work has been published recently on the combination of ultrafast imaging and ST, for the assessment of blood velocity [28,29,38] and vessel wall strain [39,40,41,42]. Here, subsample resolution displacement estimates were obtained by performing 1D [28,29,38,43] or 2D [42] parabolic interpolation of the pattern matching function. No clear reasoning was provided concerning the match between this interpolation function and the sampling, or shape, of the pattern matching function, which is a direct result of the spacing in the ultrasound sampling grid.

In this work, we propose a new concept that combines the choice of the beamforming grid with the choice of the subsample estimator. By doing so, the available data and interpolator can be matched. For this, the sampling grid for ultrafast data was based on the dimensions of the point-spread-function (PSF) of the imaging system. When ST is performed, the pattern matching function is hereby captured in a standardized way, with equal amount of information, i.e., sample points, in the axial and lateral direction throughout the peak. The use of a 2D cubic subsample interpolator, which matches with the shape of the peak, allows for a joint estimate of subsample displacement in both directions. We hypothesize that this will increase the accuracy of subsample displacement estimation and thereby will enlarge the range of detectable displacements and velocities. The concept of this technique and first simulation results on parabolic blood flow were presented in a conference proceeding [44]. The present paper fully describes the reasoning and theory behind the method and presents in vitro rotating disk experiments and carotid artery vessel wall simulations to further study the performance of the method for both vector velocity imaging and strain estimation in the carotid artery. Both wall and flow dynamics were studied, showing a wide range of inter-frame displacements and displacement gradients, to investigate the applicability of the proposed method in a broad field.

2. Materials and Methods

2.1. Theory

2.1.1. Multi-Step Speckle Tracking

Tissue displacements were estimated using a multi-step ST algorithm. Such multi-scale methods are able to obtain accurate displacement estimates in highly discontinuous displacement fields [45,46]. Normalized 2D cross-correlation was used as the pattern matching function. Inter-frame displacements were first estimated using relatively large 2D windows of envelope (demodulated RF) data, where after, in the second step of the algorithm, smaller windows of RF data were used to refine the displacements. Subsample displacement estimates were obtained by performing interpolation of the cross-correlation-function (CCF).

2.1.2. PSF-Shape-Based Beamforming Grid

The shape of the peak of the CCF is related the PSF of the imaging system, i.e., the combination of transducer choice and transmit sequence. Together with the data sampling grid, the PSF strongly influences the appearance of the peak in the CCF. By incorporating knowledge about the PSF dimensions into the data sampling grid, the CCFs can be captured in a standardized way, with an equal number of samples in the axial and lateral direction throughout the peak. These so-called circularly shaped peaks, i.e., peaks with circular isocontours, provide an equal amount of sample points to the subsample estimator in both directions. This ensures the interpolator gives equal weight to the available information in both directions. Using a 2D interpolator, which matches with these circular shaped peaks, allows for joint 2D subsample displacement estimation. In this work, we investigate whether standardization of the shape, or sampling, of the peak combined with 2D cubic interpolation increases the accuracy and precision of multi-step ST-based displacement estimation.

Circular-shaped CCF peaks were realized by choosing the data sampling grid such that a single period of the system’s PSF contained an equal amount of sampling points in axial and lateral direction. Firstly, concerning the RF data (Figure 1a), six samples per wavelength were beamformed in axial direction, in accordance with the all-positive-samples condition [34], resulting in an axial sampling distance of d_ax,rf = c/(2 × f_c × 6), with c the speed of sound and f_c the center frequency of the emitted waveform. The lateral sampling distance was determined according to d_lat,rf = d_ax,rf × PSF_ratio × N_cycles, with PSF_ratio the average ratio between the lateral and axial PSF dimensions and N_cycles the number of cycles present in the PSF axially. The envelope data sampling grid (Figure 1b) was generated by down sampling the RF sampling grid in axial direction with the factor N_cycles, resulting in an axial spacing of d_ax,env = d_ax,rf × N_cycles. The lateral spacing for the envelope data was set similar to the lateral spacing of the RF data (d_lat,env = d_lat,rf). With these spacing settings, the peak of the CCF will appear circular-shaped when using both envelope and RF data in the multi-step ST algorithm.

2.2. Imaging Setup

An L5-12 ATL (ATL, Bothell, WA, USA), 38 mm linear array transducer was used in this work. Zero-degree plane waves were transmitted by simultaneously activating 128 elements at the center of the transducer, using a rectangular apodization window. The same 128 elements were also active in receive. The transducer parameters are listed in

Table 1. Transducer properties.

Properties L12-5, 38 mm	Experiments	Simulations
Center frequency (fc)	8.9 MHz	9 MHz
Sampling frequency	36 MHz	180 MHz
Transducer element pitch	198 µm
Transducer element height	4 mm
Transducer element width	173 µm
Nr. of transducers elements	192
Nr. of active elements	128
Excitation signal	3-cycle pulse at fc

. The Field II ultrasound simulation program [47,48] was used for simulations and the Verasonics Vantage 256 experimental ultrasound system (Verasonics, Kirkland, WA, USA) for the experimental measurements.

2.3. Simulations

Ultrasound RF element data of blood and tissue were simulated, which are represented by a set of random point scatterers. About 300 scatterers per mm³ where used and each physical transducer element was subdivided into 10 by 10 mathematical elements to increase computational accuracy. A sampling frequency of 180 MHz was used, the resulting RF element data were down sampled to 36 MHz to mimic experimental settings. The speed of sound was assumed to be 1540 m/s.

2.3.1. Parabolic Flow in Vessel Phantom

Parabolic flow through a rigid, straight tube with a peak velocity (v_max) of 1.5 m/s was simulated at flow angles of 75° till 105°, with steps of 3°. The vessel had a lumen diameter of 6 mm and was centered at 20 mm depth. Scatterers, representing the blood, had Gaussian distributed amplitudes (mean = 5, SD = 1). Plane wave ultrasound data were simulated for pulse repetition frequencies (PRFs) of 2.5 till 15 kHz, with steps of 2.5 kHz.

Table 2. Imaging pulse repetition frequency (PRF) and resulting inter-frame (IF) displacements.

Simulations	PRF	IF Displacement Range (µm)
Blood flow	15 kHz	0–100
	2.5 kHz	0–600
Vessel wall	25 Hz	0.5–11.5
	50 Hz	2.5–72
Rotating Disk Experiments	PRF	IF Displacement Range (µm)
Blood flow scenario	8 kHz	0–23
	4 kHz	0–46
	2 kHz	0–91
	1 kHz	0–182
Vessel wall scenario	8 kHz	0–23
	4 kHz	0–46

provides an overview of the simulated PRFs and resulting inter-frame displacements. For each unique combination of flow angle and PRF, a pre and post frame were simulated, where in between the scatterers were propagated according to the parabolic velocity profile. This process was repeated ten times, with random initial scatterer positions, allowing 10 independent velocity estimates. No clutter filtering or time averaging were performed. Band limited noise was added to the beamformed RF data resulting in a signal-to-noise ratio (SNR) of 12 dB.

2.3.2. Carotid Vessel Wall Displacements

Simulations of the carotid artery vessel wall were performed using deformations patterns as obtained from a patient-specific finite element model (FEM) of the carotid artery at the bifurcation [40]. A transversal image plane, containing the internal and external carotid artery, was simulated. The amplitude of the scatterers representing the carotid artery were set at two, the surrounding tissue scatterers had an amplitude of one. Specular reflections were mimicked by positioning scatterers at the lumen-vessel wall and vessel wall-surrounding tissue transitions, which enhanced the realism of the RF data. Simulations of the vessel wall were performed in early systole, where peak inter-frame displacements are present and in late diastole, where low inter-frame displacements are present. A pre and post deformation frame were simulated for each time point, using a PRF of 50 Hz and 25 Hz (see Table 2) in systole and diastole, respectively. Band limited noise was added to the beamformed RF data resulting in an SNR of 30 dB.

2.4. Rotating Disk Experiments

A rotating disk experimental setup was built, providing a range of velocities at all possible beam-to-flow angles, or, a range of inter-frame displacements at full 360° range of directions. This type of setup has been used to study the performance of velocity estimation methods [29,43].

A cylindrical homogeneous phantom was created using a 15% polyvinyl alcohol (PVA) solution [39]. The solution was poured into a cylindrical mold (Øinner = 20 mm), containing a fixating-structure (see Figure 2) to attach the phantom onto the motor unit during the experiments, and subjected to three freeze-thaw cycles. Shrinkage of the phantom was present during assemblage, resulting in a final diameter of 18.2 mm. The phantom was attached to a motor unit (Closed Loop Step Motor, ARM69AC, Oriental Motor USA Corp., Torrance, CA, USA) and placed in a water tank (see Figure 2). It was rotated at a constant angular velocity, resulting in a maximum velocity at the outer edge of 18.2 cm/s. The phantom was positioned perpendicular to the transducer’s scan plane at a depth of 3.4 cm. Plane wave data were acquired continuously using a PRF of 8 kHz for 4 s.

Two different approaches were used for processing the acquired data. One approach resembling flow velocity estimation and the second approach analogous to the processing performed for the vessel wall displacement simulations.

2.4.1. Rotational Flow

After RF data beamforming and prior to velocity estimation, clutter filtering was performed. Although the phantom does not contain actual blood or tissue, including clutter filtering will give a more realistic picture of the performance of the methods for velocity estimation. A finite impulse response (FIR) filter with an order of 46 was used in combination with a temporal sliding window, not sacrificing frames for filter initialization. The −3 dB velocity cut-off point was equal to 0.068 times the Nyquist velocity. By down sampling the beamformed RF data in temporal direction, a range of PRFs (8, 4, 2 and 1 kHz) was studied. The resulting inter-frame displacements are summarized in Table 2. The same filter characteristics were used for all PRFs and estimated velocities were median averaged over 40 temporal frames.

2.4.2. Rotational Vessel Wall Displacement

Inter-frame displacements were estimated after RF data beamforming. Both PRFs of 8 and 4 kHz were studied, resembling inter-frame displacements present in early systole and late diastole of the vessel wall dynamics (see Table 2).

2.5. Beamforming Setup and Grid Definition

To be able to calculate the spacing for the PSF-shape-based beamforming grid, the PSF dimensions were determined for both the simulation and experimental setup (see

Table 3. Beamforming parameters. PSF: point-spread-function.

Grid Type	Parameter	Simulations	Experiments
PSF-shape-based grid	PSFratio	1.5	2.0
	dax,RF	14.2 µm	14.4 µm
	dlat,RF & dlat,ENV	107 µm	143 µm
	dax,ENV	71 µm	72 µm
Conventional grid	dax	21.4 µm	21.6 µm
	dlat	198 µm	198 µm
Both grids	Apodization window	Hamming
	F-number	0.875

). Field II was used to simulate the response of multiple point scatterers placed at spatial positions relevant for carotid artery imaging (Figure 3a). Simulated element data were beamformed on a high resolution grid and intensity curves for each individual scatterer were generated (Figure 3b). The axial (PSF_ax) and lateral (PSF_lat) size of the PSFs were determined based on the −6 dB values. The PSF_ratio was determined according to PSF_ratio = PSF_lat/PSF_ax. The resulting ratios were averaged for all simulated scatterers. The PSF dimensions were determined experimentally using a wire phantom (60 μm diameter) placed in a water tank. Plane wave data were acquired of the wire phantom and similar processing was performed to determine the PSF_ratio.

Delay-and-sum beamforming was used to generate RF data for every plane wave transmission at pre-specified data sampling points [39]. The data were beamformed on the PSF-shape-based grid using the axial and lateral sampling distance as specified (Section 2.1.2) for the RF data (d_ax,RF and d_lat,RF). Consequently, the envelope data were obtained by demodulation and down sampling of the RF data in the axial direction. Furthermore, the element data were also beamformed on a conventionally spaced grid, with a lateral sampling distance equal to the pitch and 4 samples per wavelength in axial direction. Details on the grid spacing and beamforming parameters can be found in Table 3.

2.6. Multi-Step Speckle Tracking Settings

Inter-frame displacements were estimated using data beamformed on both the conventional grid and the PSF-shape-based grid. The size of the search windows was tuned to the maximum displacement expected in each situation. Spatial smoothing of the estimated displacements was performed using median filters. Inter-frame displacement estimates were multiplied by frame rate to calculate velocities.

Subsample resolution was resolved by interpolation of the CCF in both steps of the algorithm. Two types of interpolation were used: 1D parabolic fitting in each direction of the CCF independently, since it is one of the most often used interpolators in ultrasound applications as described in the Introduction, and 2D cubic interpolation, where a joint estimation of subsample displacement in both directions was obtained. The parabolic fitting is solved analytically using the commonly applied three-point parabola fitting [25,33]. 2D cubic interpolation was performed with a resulting effective upsampling factor of 10 and 1000 in the first and second step of the algorithm, respectively. For the second step, a fast in-house built implementation was used, where the high final resolution was reached iteratively by zooming in onto the peak of the CCF in multiple iteration steps. Given the CCF with its maximum at position (i,j), the new sampling points (ix,iy) are defined by ix = ixm*5^iter + i and iy = iym*5^iter + j, with matrices ixm and iym representing the direct neighbors of (i,j) according to ixm = [−1 0 1;−1 0 1;−1 0 1] and iym = [−1 −1 −1; 0 0 0; 1 1 1]. Ten iteration steps (iter = 1–10) were used, resulting in the effective upsampling factor of 1000. A complete overview of the settings used in the ST algorithm can be found in

Table 4. Settings for 2-step speckle tracking algorithm.

Parameter	Step	Blood Flow Simulations	Vessel Wall Simulations	Blood Flow Experiments	Vessel Wall Experiments
Template1	1	0.84 × 2.20	0.93 × 0.77	0.9 × 1.6
-	2	0.63 × 2.00	0.32 × 0.38	0.42 × 1.0
Search window1	1	2.31 × 8.14	1.01 × 1.63	1.34 × 3.0
-	2	0.74 × 2.85	0.43 × 1.24	0.54 × 2.14
2D median filter2	1 & 2	5 × 3 3	5 × 3 3	12 × 8 3
-	-	5 × 5 4	5 × 5 4	12 × 11 4

¹ Axial × lateral window size (mm); ² Axial × lateral window size (samples); ³ Settings used for data on conventional grid; ⁴ Settings used for data on PSF-shape-based grid.

.

Summarized, four different combinations of beamforming grid and CCF interpolation type were studied in this work: (1) data beamformed on the conventional grid combined with 1D parabolic interpolation (conv_1Dpar) and (2) combined with 2D cubic interpolation (conv_2Dcub), (3) data beamformed on the PSF-shape-based grid combined with 1D parabolic interpolation (PSF_1Dpar) and (4) combined with 2D cubic interpolation (PSF_2Dcub).

The performance of the four approaches was studied in terms of the accuracy and precision of the velocity magnitude, velocity angle and the displacement estimates. The accuracy of the velocity magnitude and the inter-frame displacement estimates was assessed by calculating the absolute error percentage between the mean estimates and the ground truth as follows:

$$\mathrm{Error} \left(\mathrm{i}\right) =\frac{\left|ES{T}_{mean}\left(i\right)-GT\left(i\right)\right|}{GT\left(i\right)}\ast 100\%$$

(1)

with i all estimated samples within the region-of-interest (ROI), EST_mean the mean of the estimates and GT the ground truth value. EST_mean was based on n = 10 realizations for the simulated parabolic flow, n = 100 ensemble averaged velocity estimates for the experimental rotational flow, and n = 100 inter-frame displacement estimates for the rotational displacement experiments. The error for the velocity angle estimates was calculated as the absolute error, i.e., error (i) =|EST_mean(i) − GT(i)|.

The precision for the velocity magnitude and inter-frame displacement estimates was assessed by calculating the relative standard deviation:

$$\mathrm{SD} \left(\mathrm{i}\right) =\frac{{\sigma }_{EST}\left(i\right)}{G{T}_{max}}\ast 100\%$$

(2)

with σ_EST the standard deviation of the same n estimates, and GT_max the maximal ground truth value. The precision of the velocity angle estimates was calculated as the standard deviation (σ_EST). For the analysis of the experiments, the center of the disk (r < 1 mm) was excluded from the calculations to avoid velocities or displacements close to zero.

The performance for all methods for the vessel wall simulation was quantified by calculating the root-mean-squared error (RMSE) between estimated and ground truth axial and lateral displacements:

$$\mathrm{RMSE} =\sqrt{\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left(EST\left(i\right)-GT\left(i\right)\right)}^2}$$

(3)

with N the number of estimated samples within the ROI, EST(i) the estimated value for sample i, and GT(i) the ground truth value for sample i.

3. Results

3.1. Simulations

3.1.1. Parabolic Flow in Vessel Phantom

In Figure 4, the absolute bias and SD values are visualized as a function of the PRF for the parabolic flow simulations. The statistics were calculated while taking into account all simulated beam-to-flow angles. When the PRF increases, i.e., when the inter-frame displacements become smaller (see Table 2), the bias and SD for both the magnitude and angle increase. This effect is most dominant for the conventional grid, while the PSF-shape-based grid shows a more stable performance over the PRF range. Both PSF_1Dpar and PSF_2Dcub show high accuracy and precision, i.e., low bias and SD values, for the entire range of PRFs. For both methods, the median magnitude bias is below 10% and the angle bias remains below 0.55°.

Figure 5 visualizes the dependency of the performance of the methods on beam-to-flow angle. The results were obtained using a PRF of 12.5 kHz. The accuracy and precision of the velocity angle estimates show a dependency on the beam-to-flow angle, whereas the magnitude estimate is less influenced by this. When the velocity field contains larger axial velocity components, i.e., for beam-to-flow angles different from 90°, the angle bias and SD values increase, especially for conv_1Dpar and conv_2Dcub. The PSF-shape-based grid combined with either 1D parabolic or 2D cubic interpolation shows a stable performance over the range of angles, with almost no increase of angle error and SD for angles other than 90°.

Figure 4 and Figure 5 show a competitive performance for PSF_1Dpar and PSF_2Dcub. Both methods show a stable performance for all PRFs, i.e., for 0 to 600 µm inter-frame displacements (see Table 2), and for all beam-to-flow angles. However, significantly lower median SD values for both velocity magnitude and angle can be appreciated for the PSF_2Dcub method (Wilcoxon signed-rank, p < 0.05). The median magnitude SD is at maximum 5.8% for PSF_1Dpar and 4.2% for PSF_2Dcub. Furthermore, the median angle SD is maximally 0.9° and 0.6° for PSF_1Dpar and PSF_2Dcub, respectively.

3.1.2. Carotid Vessel Wall Displacements

Figure 6 visualizes the estimated inter-frame displacements in the external and internal carotid vessel wall using the four different methods. Major performance difference can be seen for the lateral displacement estimation, with PSF_2Dcub showing the highest similarity with the ground truth. The RMSE values are visualized in Figure 7. The accuracy of the displacement estimates in the axial direction is higher than in the lateral direction, showing lower RMSE values for all four methods. This can also be appreciated from Figure 6. The largest performance differences between the four methods are observed for the lateral estimates. The PSF_2Dcub method shows lowest RMSEs for both the diastolic and systolic phase.

3.2. Rotating Disk Experiments

3.2.1. Rotational Flow

Figure 8 shows the estimated and reference velocity magnitude and angle for the rotating disk experiment. Single ensemble-averaged velocity estimates for a PRF of 4 kHz are shown. Highest resemblance between velocity magnitude and ground truth is found for the PSF_2Dcub method. Here, largest deviations from the ground truth are observed in the center vertical slice of the disk, where velocities are mainly in lateral direction. This is a combined effect of the lower lateral resolution and clutter filtering. Due to the smaller axial velocity component in these regions, the estimates will be more influenced by the clutter filter [49,50]. In the bottom row of Figure 8, the angle estimation performance is visualized. Both the conv_2Dcub and PSF_2Dcub methods show good resemblance with the ground truth velocity angle.

In Figure 9, the absolute bias and SD values are visualized as a function of the PRF. Similar effects can be seen as in the parabolic flow simulations: when inter-frame displacements become smaller, the bias and SD for both velocity magnitude and angle increase. This is most profound for the conv_1Dpar and PSF_1Dpar methods. The conv_2Dcub method shows a stable performance over the range of PRFs, although slightly higher bias and SD values are found as compared to the PSF_2Dcub method, especially for a PRF of 8 and 4 kHz. The PSF_2Dcub method shows a stable performance over the PRF-range, with low bias and SD values. Even for very small inter-frame displacements, the accuracy and precision of this method are excellent. A median velocity magnitude bias of maximally 10% and a median angle bias of 9.2° are found. The median SD for the velocity magnitude and angle are maximally 4.6% and 6.5°, respectively.

3.2.2. Rotational Vessel Wall Displacements

In Figure 10, results for the rotating disk experiments are visualized when processed as vessel wall. Similar as in the vessel wall simulation study, the performance of the four methods is comparable for the axial displacement estimation. Furthermore, the axial accuracy is higher as the lateral displacement estimation accuracy. The lowest median lateral displacement error can be appreciated for the PSF_2Dcub method, although interquartile ranges of all methods show overlap. Furthermore, the precision of the lateral displacement estimation is highest for the PSF_2Dcub and conv_2Dcub methods. Overall, the PSF_2Dcub method shows the most stable performance with respect to the PRFs.

4. Discussion

In this paper, a new method has been proposed to determine the beamforming grid for ultrafast data based on the dimensions of the imaging system’s PSF. As a result of this data sampling, the CCFs used in the ST algorithm will contain circular-shaped peaks. It was hypothesized that once combined with 2D cubic interpolation, this would increase the accuracy and precision of subsample displacement estimation and thereby enlarge the range of detectable displacements and velocities. Both simulations and experiments were conducted to characterize the performance of the method. A comparison was made with conventionally spaced beamformed data and the very commonly used 1D parabolic subsample displacement estimation method. The results confirmed our hypothesis: the PSF-shape-based beamforming grid combined with 2D cubic interpolation (PSF_2Dcub) showed the most accurate and stable performance with respect to the range of inter-frame displacements, both for the assessment of blood flow and vessel wall dynamics. Besides, the performance of the method was least affected by motion direction or beam-to-flow angle.

An important finding was that largest performance differences between the four methods were found for the lateral displacement component. This can be appreciated for the vessel wall (Figure 6, Figure 7 and Figure 10), but also for the flow results (Figure 8), in the upper and lower part of the disk, where velocities are predominantly in lateral direction. Furthermore, the accuracy and precision of the axial displacement estimates is consistently higher as compared to the lateral estimates, due to the higher resolution and presence of phase information. This also explains why the performance difference is less obvious in the axial direction and very clearly visible in the lateral direction. No phase information is available and the frequency content is lower in the lateral direction. Consequently, more benefit is gained from the proposed methodology.

Several functions can be used to interpolate the CCF for subsample displacement estimation. In this work, we used 1D parabolic interpolation, one of the most widely used functions, to compare to our proposed methodology. The results however, show that 1D parabolic is less optimal for subsample displacement estimation. Large bias, SD and RMSE values were found when 1D parabolic interpolation was used for CCF interpolation. More specifically, large bias and SD values on the lateral axis were shown when subsample displacements existed on the axial axis. This effect can be appreciated in Figure 8, where vertically stripped patterns are visible in the magnitude and angle estimates for the PSF_1Dpar en conv_1Dpar methods. These patterns originate from the inaccurate lateral displacement estimates in regions of subsample axial displacements, as also reported earlier [33]. The benefit of using 1D parabolic interpolation is the low computational cost, whereas 2D cubic interpolation increases computational time. Therefore, the trade-off between computational time and performance needs to be considered.

One limitation of this study is that a single 2D interpolator was investigated. Many other studies have investigated the use of other shapes and functions, such as cosine, parabolic or spline [25,33,35,51,52], most often in 1D. However, none of these interpolators have been studied before in the context of the proposed methodology, i.e., in combination with the freedom of choosing the beamforming grid and thereby influencing the CCF shape. Further investigation should be performed to compare different 2D interpolators in combination with the PSF-shape-based grid. 2D fitting of the CCF peak instead of interpolation is not considered in this work, due to the enormous increase of computational time, which makes its use less practical.

Another method to improve subsample displacement estimation is by interpolation of the echo signal. This method was often used in combination with conventionally acquired data and reduced the error by the up sample factor. However, the computational cost increased significantly [25,53]. The benefit of the use of ultrafast data is that instead of interpolating the echo signals, actual focusing of the signals can be performed in receive, as a post-processing step. This work shows that no excessive sampling of the RF data is required when using knowledge of the PSF dimensions to beamform the data combined with a suitable 2D subsample interpolation method.

The proposed methodology is based upon the statement that the CCF shape is primarily determined by the shape of the PSF. While it is one of the most important factors influencing the appearance of the CCF, it is known that other factors will also be of influence. In general, the performance of motion estimation methods will degrade as a result of signal decorrelation. Factors such as out-of-plane motion, presence of strong gradients within the displacement field and large deformations can hamper the match of the signal template within the search window, thereby influencing the appearance of the CCF. Smartly choosing data type and window sizes is recommended in those situations [25,54] and aligning and stretching methods are suggested to minimize signal decorrelation [55,56]. Although these and other factors will also influence the CCF shape, the proposed methodology was successful under the circumstances tested in this work.

To determine the spacing of the PSF-shape-based beamforming grid, a single average PSF_ratio was used in this work. This was allowed since for carotid artery imaging, the PSF is nearly constant within the region of interest. However, for other applications, such as cardiac imaging, the size of the PSF will change with respect to the position in the imaging plane, amongst others due to the frequency dependent attenuation and location of the echo with respect to the ultrasonic beam [34]. These effects can be seen in the work of Tong et al. [57]. However, application of the PSF-shape-based grid is still possible, when using variable grid spacing tuned to the PSF shape locally. Nillesen et al. [58] showed first results using such a multi-zone, depth-dependent PSF-shape-based grid for cardiac motion imaging. The beamforming grid resolution was adjusted to the local PSF shape. This work also shows that the proposed beamforming strategy is not limited to linear array imaging. The PSF-shape-based beamforming grid can easily be used for data obtained with different transducer types, such as a phased array. The only requirement for using this method is the characterization of the imaging system by determining the (local) dimensions of the PSF.

In this work, the dimensions of the PSFs were determined using the response of single scatterers. When processing data acquired in vivo, a more sophistic approach could be favorable, which includes other factors influencing the shape of the PSF, such as attenuation. A representative estimate of the PSF_ratio could be obtained by using the spatial auto covariance method [59] to determine the speckle size in a B-mode image. Since speckle size is proportional to the dimensions of the PSF [60,61], this will provide the required information to determine the PSF-shape based beamforming grid.

In the present work, a large range of inter-frame displacements were studied by imaging with different PRFs. A stable and accurate performance of the proposed method was shown with respect to this large range of displacements. This shows the applicability of the proposed methodology for a broad field of research: it allows to study the dynamic behavior of multiple tissue types, such as the carotid artery, the heart and breast tissue, where very small displacements need to be captured, but also very large displacements and deformations can be present simultaneously.

In this work, the proposed beamforming strategy was combined with cross-correlation-based ST. However, this strategy is also applicable to other pattern matching functions, for example the sum of absolute squared differences. Furthermore, while this work presents a methodology in 2D, it could be easily adapted to 3D. The system’s PSF needs to be characterized in 3 dimensions and interpolation of the 3D CCF in 3 directions should be performed simultaneously.

Zero-degree plane waves were used to acquire data in this work. The estimation of vessel wall displacements and blood flow velocities could be further improved by performing displacement compounding [40,49,62], a method introduced by Techavipoo et al. [63] and adapted for vascular strain estimation [64] and blood flow imaging [52]. Angled plane waves are used to derive the displacement or velocity vector, using solely axial displacement estimates. However, to improve the accuracy and precision of the estimates, often 2D ST is performed to acquire these axial displacements [62,64,65], meaning the accuracy of both the lateral and axial displacement estimates influence the performance of the compounding method. Therefore, displacement compounding would also strongly benefit from the methodology as presented in this work.

5. Conclusions

In this paper, a framework for robust displacement estimation in ultrafast ultrasound data was presented. It can be used as a protocolled way to beamform ultrafast data and obtain accurate estimates of the tissue displacements. The success of this new methodology is based on the following important concepts: the matching shape of the CCF peak and the interpolation function, the equal number of sample points in axial and lateral direction throughout the peak, and the use of 2D instead of 1D interpolation functions to acquire joint estimation of subsample displacements in both directions. The proposed method defines the beamforming grid based on the system’s PSF to generate standardized circular-shaped CCFs which match with a 2D cubic interpolation function. The results presented in this work confirm the hypothesis; increased accuracy and precision were found for 2D motion estimation using the PSF-shape-based beamforming grid combined with 2D cubic subsample interpolation of the CCF. This was shown for a wide range of displacements or velocities present in carotid vessel wall and blood flow dynamics.