Shear Wave Velocity Estimation Using the Real-Time Curve Tracing Method in Ultrasound Elastography

: The estimation of shear wave velocity is very important in ultrasonic shear wave elasticity imaging (SWEI). Since the stability and accuracy of ultrasonic testing equipment have been greatly improved, in order to further improve the accuracy of shear wave velocity estimation and increase the quality of shear wave elasticity maps, we propose a novel real-time curve tracing (RTCT) technique to accurately reconstruct the motion trace of shear wave fronts. Based on the curve ﬁtting of each frame shear wave, the propagation velocity of two-dimensional shear waves can be estimated. In this paper, shear wave velocity estimation and shear wave image reconstruction are implemented for homogeneous regions and stiff spherical inclusion regions with different elasticity, respectively. The experimental result shows that the proposed shear wave velocity estimation method based on the real-time curve tracing method has advantages in accuracy and anti-noise performance. Moreover, by eliminating artifacts of shear wave videos, the velocity map acquired can restore the shape of inclusions better. The real-time curve tracing method can provide a new idea for the estimation of shear wave velocity and elastic parameters.


Introduction
The elastic information of biological tissue can provide valuable references for the diagnosis of tissue function and lesions [1][2][3][4][5]. In 1998, Sarvazyan et al. proposed ultrasonic shear wave elasticity imaging technology [6]. Sandrin et al. utilized ultrafast ultrasonic imaging to catch the resulting shear wave and proposed the real-time shear wave elastic imaging algorithm in 1999 [7], which promotes the rapid development of ultrasonic elastography. There is abundant evidence that shear wave elastography (SWE) can effectively and reliably quantify the elastic characteristics of various tissues [8][9][10][11]. Ultrasonic shear wave elastography has developed rapidly because it can calculate the elastic distribution of soft tissue according to the propagation velocity of shear waves. This technology includes transient elastography (TE) [12], acoustic radiation force impulse imaging (ARFI imaging) [13,14], comb-push ultrasound shear elastography (CUSE) [15,16], supersonic shear imaging (SSI) [17,18], etc. Supersonic shear imaging is one of the most advanced shear wave elastic imaging methods, and it shows great advantages in noninvasive evaluations of liver and breast [19][20][21]. In addition, it is widely used in clinical trials [22,23], which can better reflect the elastic information of tissue.
When using the supersonic shear imaging technology, there are a large number of techniques for reconstructing the velocity image based on the mechanical properties. At present, the elastic value of tissue is characterized by estimating shear wave velocity, including the time-to-peak (TTP) method [24], the correlation-based method [25], the lateral peak (LP) method [26], the algebraic inversion method [27], Radon transform method [28,29], etc. However, the uncertainty of noise in ultrasonic motion signal limits the aforementioned methods' application in the ultrasonic direction [26]. Based on this, McLaughlin et al. proposed a method to estimate the local shear wave velocity by obtaining the displacement information of shear waves through time signal cross-correlation in 2006 [25]. Palmeri et al. proposed the time-to-peak (TTP) method to characterize the shear wave velocity by using the time peak displacement data of the acoustic radiation force response outside the excitation region in 2008 [24]. The time-to-peak method has been successfully applied to shear wave elasticity imaging [30,31]. In order to improve the robustness of the time-to-peak method, Wang et al. introduced an application of random sample consensus (RANSAC) to estimate shear wave velocity from tissue displacement tracked by ultrasound in vivo [32]. Carrascal et al. proposed two modified shear wave velocity estimation algorithms in 2017, one of which is to search for the maximum motion in space and time, and the other of which is to apply amplitude filters to select points with a motion amplitude higher than the threshold of the shear wave group velocity estimation [26]. As this kind of peak method depends on the accurate identification of peak information, it is still necessary to improve the method of shear wave velocity estimation in the presence of noisy motion data.
We propose a real-time curve tracing (RTCT) method for shear wave trajectory, and carry out velocity estimation and two-dimensional shear wave velocity image reconstruction since the aforementioned peak method is susceptible to the local maximum when estimating the shear wave velocity. Based on the peak method, the real-time curve tracing method starts from the perspective of image processing and improves the accuracy of shear wave velocity estimation to a certain extent through connected domain processing. The accuracy of the real-time curve tracing method is proved on an ultrasound phantom by using inclusions that are softer and harder than the background material for velocity image reconstruction. At the same time, the influence of variable Gaussian noise on the accuracy of the method is studied to verify its robustness. In addition, we reproduce the velocity estimation results of two different peak measurement methods in the same phantom, and make a comparison and analysis.

Method and Principle
In this paper, shear waves are generated by ultrasonic radiation force. Firstly, the relationship between the energy of shear wave and the distance of the focus point is obtained by the log compressed B-mode data, that is, to search the spatial distribution of each time i (i = [1 : i]) on all i time points, as shown in Equation (1).
where the displacement spatiotemporal data is given by u(x, t), and the arg{·} operator is intended to obtain the indices for the points related to the maximum of the displacement profiles. The position information of each point corresponding to S i is found in the frame i th image f i (x, y) of the shear wave propagation video (where (x, y) are the pixel coordinates, the transverse pixel spacing is 0.25 λ and the longitudinal pixel spacing is 0.5 λ, 1 λ = 0.2464 mm). Taking these points as the center, a range is determined with the standard of excitation duration, and then convolution is performed in the range with a 5 × 5 Sobel operator and f i (x, y). In other words, the position of each point corresponding to C i is taken as the center, and the edge is searched i th times in a certain range. The approximate luminance difference matrices G x and G y in x and y directions are obtained, as shown in Equations (2) and (3), where * is the convolution operation.
Finally, the gradient magnitude M and direction θ of each pixel are calculated. The gradient amplitude and direction angle are obtained by polar coordinate transformation, which are expressed by Equations (4) and (5).
where M(x, y) and θ(x, y) reflect the intensity and the direction of the image edge, respectively. In order to determine the edge information more accurately, we keep the maximum value of the local gradient and suppress its non-maximum value. The gradient intensity of the current pixel M(x, y) is compared with the gradient intensity of four pixels along its positive and negative gradient directions. If the gradient intensity of the current pixel is higher than that of the other four pixels, the pixel is retained as an edge point. After that, we need to select and connect the edge points, and the specific operation is as follows: (1) selecting the appropriate threshold T, (2) marking the edge points whose gradient amplitudes are greater than T, (3) connecting the marked edge points to form a closed curve.
After i times of searching, the edge information of each search is determined, and it is connected to form a closed curve. Then, the intersection of the closed curves containing any point of C i in the i times search is taken as a connected domain. Meanwhile, to eliminate the artifacts of shear wave propagation, this paper includes a directional filter matrix. A directional filter with a minimum circumscribed rectangle is designed for the connected domain of the shear wave image, and the data outside the directional filter are set to 0. Next, Gaussian filtering is applied to the shear wave image h i (x, y) after this processing, and the filtered image p i (x, y) is obtained, as shown in Equation (6).
where K(x, y, σ) is the two-dimensional Gaussian kernel function, σ is the standard deviation and * denotes convolution operation. Finally, the connected domain of the directional filter in image p i (x, y) is regarded as the shear wave front, and the center points in each transverse direction in the connected domain of the shear wave front are connected into a smooth curve to complete the real-time curve tracing of the shear wave front.
From the perspective of image processing, the real-time curve tracing method utilizes a curve to fit the discontinuous shear wave with a certain width, which can realize the real-time visualization of shear wave propagation. Based on the above principle, we reconstruct the shear wave velocity image by using the real-time tracing curve of shear wave propagation. For the selected region of interest (ROI), the shear wave velocity V(x, y) of the pixel point (x, y) is given by Equation (7).
where ∆x is the displacement of the shear wave tracing curve relative to the previous frame, and ∆t is the interval time between each frame of the shear wave. Equation (7) can be used to calculate the propagation velocity of the shear wave in the region of interest at each pixel, so as to reconstruct the shear wave velocity image. Certainly, interpolation and filtering algorithms can make the speed reconstruction image have better resolution and contrast. Then, the reconstructed shear wave velocity image is transformed into a two-dimensional ultrasonic elastic image by Equation (8).
where E is Young's modulus, ρ is the density of the tissue and c is the shear wave propagation velocity.

Experimental Environment
The excitation and acquisition of the shear wave are realized by the Verasonics ultrasonic experimental system (America), and its model is Vantage 256. The L11-4V ultrasonic linear array probe is used in the experiment, and its center frequency is 6.25 MHz. We use an ultrafast plane wave imaging technique to image the propagation of the resulting transient plane shear waves. Other parameters of the system are shown in Table 1. Four experiments are conducted in a 040GSE ultrasound phantom with a background elasticity of 30 kPa by the supersonic shear imaging method. The selected experimental areas are shown in Figure 1, in which the red rectangle represents the region of interest for each experiment, and the dotted box in the red rectangle is the range for velocity estimation. In other words, the dotted box is considered to be the main propagation area of the right shear wave in each group of experiments, which is also our targeted imaging area.
In Figure 1a, a larger solid black circle and a smaller solid black circle on the right side of the region of interest represent the anechoic stepped cylinders. The circles with different decibels on the lower side of the region of interest represent grayscale target groups, ranging from −9 dB to +6 dB. In the region of interest of Figure 1b-d, there are inclusions with a diameter of 6 mm, which are different from the background elasticity, and their elasticity is 10 kPa, 40 kPa and 60 kPa, respectively. In these figures, the small white dots far smaller than the size of the inclusions are the target points, and the white circle without markings is a hyperechoic region of the same size as the other grayscale target groups.

Shear Wave Velocity Estimation
In this paper, the real-time curve tracing is performed for the 10th frame image of shear wave propagation in the targeted imaging area of Figure 1a. Figure 2a shows the original propagation trajectory of the shear wave. The shear wave front is obtained by the directional filter matrix, and then the central point of each transverse shear wave front is connected, as shown with the red curve of Figure 2b. Since the shear wave propagates in a homogeneous region, the tracing curve is close to a straight line with an inclination angle. Then, we do the same operation for the 30th frame image of shear wave propagation in the targeted imaging area of Figure 1d. Similarly, Figure 2c shows the original propagation trajectory of the shear wave, and the tracing curve is shown with the red curve of Figure  2d. Owing to the deformation of the shear wave front at the inclusion, the tracing curve protrudes to the right, which is consistent with the expected result. Because the shear wave front cannot cover the whole longitudinal range, we set two transverse lines at the head and tail of the shear wave front. Generally, the isotropic homogeneous region is used to analyze and evaluate the estimation method of shear wave velocity [33]. For the homogeneous region, the shear wave velocity can be quickly estimated by taking the lateral moving distance of several tracing curves in the same transverse direction as the shear wave displacement.
With the real-time curve tracing method, the tracing curves of each frame of the highframe rate shear wave propagation in the homogeneous region of Figure 1a can be obtained. Figure 3 lists the six tracing curves that are initially frame 4 and frame interval 5. Directly perceived from Figure 3, the shear wave tracing curve is an approximately inclined straight line, and the interval of each curve is roughly equal. On the one hand, these characteristics reflect the Mach cone shape of shear wave and, on the other hand, they reflect that the shear wave propagates at a uniform speed, which indicates that the targeted area is homogeneous. Since the longitudinal range of the data that can be used to calculate the wave velocity is different in each frame, the common part of these frames is selected to estimate the shear wave velocity (about 25 λ to 85 λ). A to F are the intersections of the tracing curves of the 4th, 9th, 14th, 19th, 24th and 29th frames at a depth of 55 λ, respectively. By performing linear regression processing on the ordinates of these six intersections, the propagation velocity of the shear wave in the specified transverse direction can be estimated. Of course, one can choose the intersection point of any number of frames at a certain transverse position for linear regression to estimate the propagation velocity of the shear wave in the homogeneous region. The estimation error decreases with the increase in frame number. To verify and compare the performance of the proposed method, the time-to-peak and the lateral peak method are used to estimate the shear wave velocity in the homogeneous region under the same conditions, taking a homogeneous region with a Young's modulus of 30 kPa as an example. The implementation processes of the time-to-peak method and the lateral peak method are shown in Figure 4. The time-to-peak method finds out the maximum time of shear wave displacement in each transverse position. The curves in Figure 4a describe the relationship between the energy of the shear wave and the time at different displacements (these data were normalized). Therefore, the shear wave velocity can be estimated from the known displacement and the time to peak, as shown in Figure 4c. The lateral peak method is used to select the information of the shear wave energy peak point to process. As shown in Figure 4b, the curves describe the relationship between the shear wave energy and distance away from the focus point at different times (these data were normalized). In that way, the velocity of the shear wave in a certain period of time can be estimated by the known time interval and the peak displacement, as shown in Figure 4d.

Homogeneous Region
According to the real-time curve tracing method proposed in this paper, the experimental data of a high-frame rate shear wave propagation video can be obtained in the homogeneous region of the 040GSE ultrasound phantom shown in Figure 1a, with which the shear wave velocity image reconstruction is realized. The B-ultrasound image of the homogeneous area is shown in Figure 5a. In the B-ultrasound image, the thick red line at the top corresponds to the position of the probe, the white rectangle is the selected region of interest and the four red cross marks in the longitudinal direction correspond to the focus coordinates in Table 1. The real-time curve tracing method is utilized to reconstruct the shear wave velocity image of the selected region of interest in the homogeneous region of the 040GSE ultrasound phantom. In order to facilitate the presentation of the reconstruction effect, the results in Figure 5b are selected for subsequent detailed analysis, with a transverse range of 20 λ to 60 λ and a longitudinal range of 30 λ to 80 λ.

Stiff Spherical Inclusion Regions with Varying Elastic Values
In this subsection, the shear wave velocity images of the region of interest containing a stiff spherical inclusion in Figure 1b-d are reconstructed, of which the B-ultrasound images are shown in Figure 6(a-1,a-2,a-3), respectively. The elasticity of the inclusions is 10 kPa, 40 kPa and 60 kPa, and the background elasticity is 30 kPa. When reconstructing the shear wave velocity image of the stiff spherical inclusion regions, the image processing is the same as that of the homogeneous area. Two targeted regions are intercepted for subsequent analysis, as shown in Figure 6(b-1,b-2,b-3). The background is selected as the black rectangle area with a transverse range of 30 λ to 50 λ and a longitudinal range of 25 λ to 35 λ. The inclusion is selected as the black circular area with a diameter of 24 λ, about 6 mm.

Noise Adaptability
In theory, even under different noises, the motion state of a shear wave can be accurately characterized as long as the propagation trajectory of the shear wave can be identified. In order to further analyze the noise adaptability of the real-time curve tracing method, the shear wave propagation images are processed under different noises. Four frame tracing curves with an initial frame number of 10 and frame interval of 5 are selected as raw reference data to conduct two sets of experiments. In the first experiment, these raw reference data are contaminated with zero-mean Gaussian noise, of which the standard deviation (STD) is set to 0.1, 0.2 and 0.3. In the second set of experiments, the mean value of Gaussian noise added in the raw reference data is set to 0.05, 0.1 and 0.2, and the standard deviation is 0.1. The results of the above two sets of experiments are given in in Figures 7 and 8, respectively.

Evaluation of Shear Wave Velocity
In the experiment, the real-time curve tracing method is used to estimate the transverse propagation velocity of a shear wave in the selected homogeneous region, and the results are compared with those of the time-to-peak method and the lateral peak method. According to the experimental method in Section 3.2, these three methods selected 50 points in each transverse direction for data analysis to ensure consistency and fairness. In the longitudinal range of 25 λ to 85 λ, 61 groups of transverse data were selected at intervals of 1 λ. With regard to the real-time curve tracing method, linear regression is performed on each group of transverse data, as shown in Figure 3, and then the average value of the 61 groups is obtained as the shear wave velocity in the homogeneous region. With regard to the time-to-peak method and the lateral peak method, the shear wave velocity is calculated by the peak information of the shear wave. The confidence interval is 95%, and the results are shown in Figure 9. The expression of velocity is average ± standard deviation, and the results are shown in Table 2. In this paper, the calculation formula of relative error X R is shown in Equation (9).
where V T is the theoretical velocity and V E is the estimated velocity.

Result of Velocity Reconstruction Image
In the homogeneous region, the experimental data of velocity image reconstruction in the intercepted homogeneous area shown in Figure 5b are quantitatively analyzed. The estimated velocity of the shear wave in the homogeneous area with the real-time curve tracing method is 2.9026 m/s, the standard deviation is 0.1032 m/s and the relative error is 8.21%.
In the stiff spherical inclusion regions with varying elastic values, the experimental data of velocity image reconstruction in the black rectangle intercepted area and the background area shown in Figure 6b are quantitatively analyzed. The comparison results are shown in Table 3. Table 3. Comparison of shear wave velocity estimation. Analysis of velocity image reconstruction in the intercepted area and the background area with the real-time curve tracing method.

Relative Error
(%) When the elasticity of inclusion is 40 kPa, 10 kPa and 60 kPa, the estimated relative errors are 8.25%, 7.05% and 6.77%, respectively, which means that the relative error of estimation velocity gradually decreases with the increase in the elastic difference between inclusion and background.

Noise Adaptability Results
From Figures 7 and 8, we can see that even in the case of noise, the real-time curve tracing method can still express the shear wave trajectory accurately. Since the longitudinal range of the shear wave tracing curve has changed after adding noise, this paper only selects and analyzes the common longitudinal range of four tracing curves, namely, 25 λ to 85 λ. Hereon, the real-time curve tracing method is used to estimate the shear wave velocity under the condition of no artificial noise or different noises, and the standard deviation and relative error are compared. The results are shown in Figure 10, Figure 11 and Table 4, and the theoretical shear wave velocity is 3.1623 m/s.   The results in Figures 7 and 8 state clearly that the change of the Gaussian noise standard deviation will affect the blur degree of the image, and the change of the mean value will affect the brightness and darkness of the image. Compared with the blur degree of the shear wave image, the brightness and darkness of the shear wave image can affect the shear wave velocity estimation more. It can be seen in Figures 10 and 11 and Table 4 that the mean value of Gaussian noise has a great influence on the shear wave velocity, which makes the standard deviation of the estimated velocity larger. Therefore, it is necessary to avoid the change of shear wave image brightness and darkness in the estimation of shear wave velocity.

Discussion
In the homogeneous region with an elastic modulus of 30 kPa, the estimated velocity, the standard deviation and the relative error with the time-to-peak method are 2.8078 m/s, 0.1944 m/s and 11.21%, respectively. The corresponding estimated values of the lateral peak method are 2.8622 m/s, 0.2114 m/s and 9.49%. However, the estimated velocity, the standard deviation and the relative error with the real-time curve tracing method are 2.9382 m/s, 0.1221 m/s and 7.09%, respectively. Evidently, compared to the time-to-peak method and the lateral peak method, the real-time curve tracing method improves the accuracy while estimating the shear wave velocity in the homogeneous area, which is due to the improved frame rate stability of hardware devices. The main factor affecting the accuracy of the real-time curve tracing method is the stability of the frame rate. If the frame rate stability is evaluated and modified, the accuracy of shear wave velocity estimation will be further improved. The experimental results show that the real-time curve tracing method can provide new ideas for the evaluation of shear wave velocity and elastic parameters.
It can be concluded from the experimental results of the homogeneous region of the 040GSE ultrasound phantom that the real-time curve tracing method has more advantages in image reconstruction of shear wave velocity because of its high accuracy in two-dimensional shear wave velocity estimation. In the results obtained by using the stiff spherical inclusion regions, the results fully illustrate that for the stiff spherical inclusion regions, regardless of the elasticity and location of inclusions, the shear wave velocity image can be reconstructed by using the real-time curve tracing method. The above experimental results verify the feasibility of the real-time curve tracing method for shear wave velocity image reconstruction in the stiff spherical inclusion regions.
In the noise experiment, the experimental results show that the influence of Gaussian noise on velocity estimation is small, which indicates that this method has strong noise adaptability. Since the time-to-peak method and the lateral peak method use peak data to calculate shear wave velocity, and peak value is easily affected by noise, the real-time curve tracing method has better anti-noise performance than the time-to-peak method and the lateral peak method.
In recent years, most methods are based on time domain data to estimate shear wave velocity. These methods estimate shear wave velocity by the characteristics of shear wave motion [34,35]. Among them, when the time signal cross-correlation is used to calculate the shear wave velocity [36], the size of the window will affect the measurement results of shear wave velocity. However, the proposed real-time curve tracing method avoids the influence of window size on shear wave estimation velocity by image processing technology.
At present, there are some methods to estimate shear wave velocity using frequency domain data, which are common methods proposed in recent years [37][38][39]. Most of these estimation methods have limited spatial range in the transverse and longitudinal directions. In order to change this defect, many research teams have extended these methods, such as making them generate images of phase velocity or viscoelastic characteristics [40,41]. In 2018, Kijanka and Urban proposed local phase velocity imaging to create phase velocity images at different frequencies in homogeneous and heterogeneous media using Fourier decomposition and moving windows [38]. One of the limitations of the local phase velocity imaging method is that the results may be corrupted by outliers or motion records with high noise content. However, the real-time curve tracing method can avoid outliers in image filtering.

Conclusions
In this paper, we present a new real-time curve tracing method to estimate the shear wave velocity. Quantitative velocity images have been achieved, particularly in the homogeneous region and in the stiff spherical inclusion regions. Compared with the traditional time-to-peak method and lateral peak method, the relative error and standard deviation of the real-time curve tracing method for velocity estimation in a homogeneous area are smaller. At the same time, the proposed method can reconstruct the shear wave velocity image more accurately. In the experiment of changing the standard deviation and mean value of Gaussian noise, the proposed method shows some robustness. The experimental results and analyses show that the real-time curve tracing method can trace the shear wave propagation in real time, and outperforms other methods in terms of the estimation accuracy. Therefore, the velocity image reconstruction and elasticity map inversion based on the real-time curve tracing method can provide more accurate and important elastic information for tissue characteristics.