2-D Minimum Variance Based Plane Wave Compounding with Generalized Coherence Factor in Ultrafast Ultrasound Imaging

Plane wave compounding (PWC) is an effective modality for ultrafast ultrasound imaging. It can provide higher resolution and better noise reduction than plane wave imaging (PWI). In this paper, a novel beamformer integrating the two-dimensional (2-D) minimum variance (MV) with the generalized coherence factor (GCF) is proposed to maintain the high resolution and contrast along with a high frame rate for PWC. To specify, MV beamforming is adopted in both the transmitting aperture and the receiving one. The subarray technique is therefore upgraded into the sub-matrix division. Then, the output of each submatrix is used to adaptively compute the GCF using a 2-D fast Fourier transform (FFT). After the 2-D MV beamforming and the 2-D GCF weighting, the final output can be obtained. Results of simulations, phantom experiments, and in vivo studies confirm the advantages of the proposed method. Compared with the delay-and-sum (DAS) beamformer, the full width at half maximum (FWHM) is 90% smaller and the contrast ratio (CR) improvement is 154% in simulations. The over-suppression of desired signals, which is a typical drawback of the coherence factor (CF), can be effectively avoided. The robustness against sound velocity errors is also enhanced.


Introduction
Owing to the high efficiency, the low cost, and the non-invasive characteristics, ultrasound imaging has been an adequate technique for medical diagnosis [1][2][3]. However, the commonly adopted line scan mode for ultrasound imaging is limited to a low frame rate since it requires a series of focused beams for a single emission, which makes it difficult for the visualization of rapid tissue motions [4]. To realize ultrafast ultrasound imaging, plane wave imaging (PWI) is proposed, which utilizes a single pulse generation rather than a series of focused beams [5]. Due to the lack of focusing on pulse emissions, the frame rate is significantly increased while the imaging quality suffers a great degradation [5][6][7]. In addition, the unfocused beams cannot offer sufficient excitation energy for the region of interest (ROI). As a result, the signal-to-noise ratio (SNR) for the backscattered echo is lower than conventional methods.
To handle the shortcomings of PWI, plane wave compounding (PWC) is proposed to achieve a compromise between the high frame rate and the satisfactory imaging quality. The original concept can date back to 1981 [8]. Then, Lu and Cheng first introduced the spatial compounding of several signal of each array element. As a result, the correlation and coherence information in the 2-D dataset is thoroughly utilized.
The remainder of this paper is divided into several sections. Section 2 introduces the background of existing techniques. Section 3 illustrates the algorithm of the proposed method in detail. Section 4 presents the simulated and experimental setup, along with the results of the proposed method and conventional methods. Section 5 provides the discussion and the final conclusion is given in Section 6.

Mathematic Model of the Plane Wave Compounding
The PWC recombines a series of plane waves with different steering angles. The received echo signals can be recorded into a 2-D data matrix. One dimension stands for different array elements, which is called the receiving aperture. Another stands for different emissions, which can be regarded as the transmitting aperture [28]. If an M-element transducer array is used for ultrasound transmitting and receiving while N steered plane waves are fired in total, an N × M size echo data matrix could be obtained. After introducing an appropriate time delay into each channel, the final matrix can be expressed as follows: where x i,j (n) is the echo signals recorded by the jth element when the ith steered plane wave is used for transmitting. n is the time step. In the conventional delay-and-sum (DAS) algorithm [29,30], all channel data is averaged through rows and columns for each imaging point: where ω R is the receiving apodization window and ω T is the transmitting one.

Minimum Variance Beamformers
The MV beamformer aims to minimize the output energy subject to the desired signal being undistorted [17]. For each single emission, the MV equation could be written as follows: where w(n) is the MV weighting vector and (·) H denotes the conjugate transpose. While x(n) = [x 1 (n), x 2 (n), . . . , x M (n)] T ((·) T is the matrix transpose) represents the signals recorded by different array elements with appropriate time delays. The MV process can be expressed as an optimization problem: since the x(n) is appropriately delayed, the steering vector, d, is equal to an all one vector. Here, R(n) is the covariance matrix. The solution to (4) is: The covariance matrix, R(n), is usually unknown and must be estimated: To enhance the accuracy of the covariance matrix estimation, the spatial smoothing technique is proposed [31]. The transducer array is divided into several overlapped subarrays. Then, the covariance matrix of each subarray is averaged to decorrelate the highly correlated array signals: where x l = [x l (n), x l+1 (n), . . . , x l+L−1 (n)] T represents the lth subarray with a length of L. According to Refs. [31][32][33], the subarray length, L, should be smaller than M/2, a half of the array size. In addition, to avoid a singular covariance matrix, diagonal loading is also necessary. It is suggested to introduce a constant into the matrix: R DL (n) = R(n) + ε·I, where I is the identity matrix and the constant, ε, is often set to be ∆ (smaller than 0.1) times of trace, R [16,17]. The modified weighting vector, w MV , can be calculated by Equation (5) using the R DL (n). Then, the beamformed output should also be averaged over subarrays, which means that Equation (3) should be rewritten into: To further reduce the noise and increase the SNR, the ESBMV beamformer is proposed [18]. The covariance matrix, R, is eigen-decomposed to project the w MV onto the signal subspace. This process can be illustrated as follows: where Λ = diag[λ 1 , λ 2 , . . . , λ K ] represents all eigenvalues (λ 1 ≥ λ 2 ≥ . . . ≥ λ K ) and v k is the eigenvector corresponding to λ k . Eigenvectors corresponding to first largest eigenvalues, which are larger than λ thre , are selected to build the signal subspace: The ESBMV weight is calculated as: Then, the ESBMV output can be acquired using Equation (8).

Coherence Factor Beamformers
The CF is defined as the ratio between the coherent sum and the incoherent sum: The GCF is a more efficient approach for the spatial spectrum based beamforming [23]. First, the M-point discrete Fourier transform (DFT) is carried out to obtain the Fourier spectrum. Usually, this process is implemented using the fast Fourier transform (FFT): where k = 0 to M − 1 is the spatial frequency index. Then, the GCF can be defined as the ratio of the low-frequency energy to the total energy: GCF = energy within a low − f requency region total energy Usually, the low-frequency region is selected by a cutoff frequency, M 0 (from −M 0 to M 0 ). For point targets, the M 0 is suggested to be 0, which means that only the direct current component is used for the GCF. For complicated situations, M 0 should be set to 1-3 to avoid artifacts [23].

Joint Transmitting-Receiving Beamforming
The JTR beamformer is a modified version of the MV beamformer for the PWC [21]. Conventional MV methods need to calculate one weighting vector for each firing event, thus the computational load would be extremely high. The JTR only needs two weights in total, the transmitting one and the receiving one, which brings a lower computational complexity (CC).
Taking account of the echo data matrix in Equation (1), the receiving covariance matrix is calculated first. Define x Rx i (n) = [x i,1 (n), x i,2 (n), . . . , x i,M (n)] T as the rows of the original matrix, which represents each single transmission. The receiving covariance matrix is estimated by averaging through not only the receiving subarrays, but also all transmitting events, which means: where T is the lth subarray of the ith transmission with a length of L 1 . The MV based weighting vector for the receiving aperture, T , can be calculated using (5) or (11).
In consideration of the ultrasound reciprocity, the transmitting aperture weighting vector can be calculated using the same procedure. Define x Tx j (n) = x 1,j (n), x 2,j (n), . . . , x N,j (n) T as the columns of the original matrix, which represents each receiving element. The transmitting covariance matrix is estimated by averaging through the transmitting subarrays and all receiving elements, as well: where x Tx j,l = x l,j (n), x l+1,j (n), . . . , x l+L 2 +1,j (n) T is the lth subarray of the array data recorded by the jth element with a length of L 2 . The weighting vector for the transmitting aperture, w Tx = [w 1 Tx , w 2 Tx , . . . , w L 2 Tx ] T , can also be obtained. To accord with the spatial smoothing technique in both the transmitting and receiving dimension, the subarray division is upgraded into the submatrix division. The submatrix with a size of L 2 × L 1 can be described as follows: where i ∈ [1, 2, . . . , N − L 2 + 1], j ∈ [1, 2, . . . , M − L 1 + 1] are the indexes of the submatrix. The JTR beamformed output of a single submatrix can be calculated by: In original JTR implementation, the final output is obtained by the averaged superposition of all z i,j (n) [21], which still leaves room for improvements.

2-D generalized Coherence Factor Weighting
Taking account of the JTR output matrix: Here, we intend to use the GCF to constrain the final output. A simple solution is to arrange all elements in Equation (19) in a row and then implement the FFT. However, since the echo data in Equation (19) is from different firing events, the 1-D FFT cannot obtain the accurate low-frequency component. Another alternative is to calculate one GCF for each firing event, which means one GCF for each row in Equation (1) [24]. This process is more accurate, but could bring a high computational load.
In this paper, we propose to use the 2-D FFT of the output matrix (19). We can acquire the 2-D spatial spectrum of the original signal: where p i,j stands for the corresponding frequency domain component. Similar to the definition of the 1-D GCF, the 2-D GCF is also calculated as the ratio between the low-frequency energy to the total energy, but this time the low-frequency component is selected in both dimensions, which means: where M 1 and M 2 are the cut-off frequency in the transmitting dimension and the receiving dimension to define the 2-D low-frequency region. Using the 2-D GCF, the final output of the proposed method could be obtained: Additional remarks should be specified here that there are two distinct differences between the modified GCF and the conventional one. First, the GCF 2 . adopts the 2-D FFT, and the 2-D low-frequency components are selected, which accords with the 2-D echo data matrix of the PWC. Second, the JTR beamformed output of all submatrices (19), instead of the original echo data matrix (1), is used as the time-domain signal of the FFT. Since different rows in Equation (1) come from different emissions, adopting a 2-D FFT on (1) could bring spatial frequency errors. Therefore, the spatial smoothing technique is adopted to divide the original matrix (1) into several overlapped sub-matrices and then these sub-matrices are MV beamformed. As a result, the new FFT on Equation (19) could bring more accuracy in the frequency domain than the conventional one, which may result in a better imaging performance.

Implementation Summary
A brief implementation summary of the proposed algorithm is given here to illustrate the method in a clear way: (1) Calculate time delays for each channel of different plane wave firing events. The echo data matrix (1) can be obtained after the calculation. (2) Compute the JTR weighting vector, w Tx and w Rx , using Equations (5), (15), and (16). The JTR beamformed output matrix (19) can be acquired from Equation (18). (3) Use the 2-D FFT to get the spatial spectrum and then calculate the 2-D GCF using Equation (21).
The final output of the proposed algorithm is obtained from (22). (4) Repeat the procedure I to III over each imaging pixel to generate the final B-mode image.

Experimental Setup
The proposed method was evaluated through simulations, phantom experiments, and in vivo studies. The Matlab simulation tool, Field II, was used to acquire simulated data [34,35]. The Verasonics ultrasound platform (V1, Verasonics, Redmond, WA, USA), as shown in Figure 1a, was used to acquire the phantom and in vivo data. For both simulations and experiments, a 5-MHz, 128-element transducer array with a 0.3 mm pitch was used (L11-4v, Verasonics, Redmond, WA, USA). The excitation pulse was a two-cycle sinusoid and the sampling rate was 40 MHz. For the PWC implementation, 49 steered plane waves at a 0.5 • interval ranging from −12 • to 12 • were emitted in total.  For B-mode image evaluations, the simulated and experimental results using the proposed JTR-MV GCF and JTR-ESBMV GCF are presented along with results of the DAS, the CF, the JTR MV, and the JTR ESBMV. First, the DAS result is shown as the reference of other methods. Then, the CF result shows the effect of only using CF-based methods while the JTR results show the influence of only using MV-based methods. At last, the results of the proposed method are given to validate the effectiveness of the integration between JTR and GCF. An additional experiment is also implemented to investigate the influence of sound velocity errors on different beamformers. A single point is simulated with an overestimation of the sound velocity by 0%, 2%, 5%, and 10%. Results of the DAS, the CF, the JTR-MV, the JTR-ESBMV, the JTR-MV GCF, and the JTR-ESBMV GCF are shown together. purpose-multi-tissue-ultrasound-phantom/. A 28-year-old male volunteer helped us to obtain the in vivo data of two human carotid arteries, a human thyroid and a human parotid.

Parameters and Evaluation Metrics
For B-mode image evaluations, the simulated and experimental results using the proposed JTR-MV GCF and JTR-ESBMV GCF are presented along with results of the DAS, the CF, the JTR MV, and the JTR ESBMV. First, the DAS result is shown as the reference of other methods. Then, the CF result shows the effect of only using CF-based methods while the JTR results show the influence of only using MV-based methods. At last, the results of the proposed method are given to validate the effectiveness of the integration between JTR and GCF. An additional experiment is also implemented to investigate the influence of sound velocity errors on different beamformers. A single point is simulated with an overestimation of the sound velocity by 0%, 2%, 5%, and 10%. Results of the DAS, the CF, the JTR-MV, the JTR-ESBMV, the JTR-MV GCF, and the JTR-ESBMV GCF are shown together.

Parameters and Evaluation Metrics
The parameters used in simulations and experiments are carefully selected. The subarray length used for the spatial smoothing is set to 0.3 times the aperture size to make a compromise between the imaging quality and the calculation time [36]. To specify, The ESBMV threshold, β, is set to 0.02 to define the signal subspace in which eigenvalues are larger than β times the largest eigenvalue [21,36]. The cut-off frequencies of the GCF are M 1 = M 2 = 0 for point targets and M 1 = M 2 = 1 for cysts and in vivo tissues according to Ref. [23]. A large cut-off frequency could avoid dark artifacts in the speckle region. The f-number is 1.4 with a rectangular window and the diagonal loading factor ∆ is set to 0.01.
For quantitative assessment, the evaluation metrics used in the following section are the full width at half maximum (FWHM) (the −6 dB bandwidth of the mainlobe) for point targets and the contrast ratio (CR), and the contrast-to-noise ratio (CNR) for cyst targets. The CR and CNR can be calculated by: where µ b and µ c are the mean magnitudes of the signal inside the speckle and cyst, and σ b and σ c are the standard deviations of the magnitude in the speckle and cyst, respectively. All transducer and acquisition parameters along with processing parameters are given in Table 1.

Simulated Study
Simulated point target images are shown in Figure 2. Figure 2a is the DAS result with visually obvious sidelobes. In Figure 2b, the CF shows its effectiveness in suppressing sidelobe noises. Simulated point target images are shown in Figure 2. Figure 2a is the DAS result with visually obvious sidelobes. In Figure 2b, the CF shows its effectiveness in suppressing sidelobe noises. Figure  2c,e are results of the MV based adaptive beamformers, which indicate that the MV can enhance the imaging resolution while the ESBMV could bring better noise reduction than the conventional MV. The JTR-MV GCF and JTR-ESBMV GCF results are presented in Figure 2d,f respectively. From the comparison between Figure 2c,e and 2d,f, it can be seen that the proposed method achieves a narrower mainlobe width than the JTR method. The narrow mainlobe and low sidelobes indicate that the proposed method could bring a high resolution. For further comparisons, lateral variation plots at z = 30 mm of different beamformers are shown in Figure 3 while the corresponding FWHMs are given in Table 2. The proposed JTR-ESBMV GCF achieves the best performance with a 90% smaller FWHM than the conventional DAS beamformer, which validates the resolution improvement by the proposed method. For further comparisons, lateral variation plots at z = 30 mm of different beamformers are shown in Figure 3 while the corresponding FWHMs are given in Table 2. The proposed JTR-ESBMV GCF achieves the best performance with a 90% smaller FWHM than the conventional DAS beamformer, which validates the resolution improvement by the proposed method.     Simulated anechoic cyst images are shown in Figure 4. In Figure 4a, noises are still at a high level inside cysts, which results from the sidelobes of background speckles. The CF result in Figure 4b obtains a satisfactory noise reduction inside cysts. However, signals inside the speckle region are also over-suppressed and the speckle patterns are severely damaged. It is due to the distortion characteristic of the CF method. The JTR method can obtain a better speckle performance and help remove the noise slightly as shown in Figure 4c,e. In comparison, the proposed method can significantly remove the noise inside cysts, which brings a more distinct cyst boundary. While the speckle patterns are still well preserved with less dark spots thanks to the modified GCF. The results confirm that the proposed method has the potential of maintaining the high contrast and good speckle performance at the same time.

Experimental Phantom Study
The wire target phantom results are shown in Figure 5, which are similar to the simulated results presented before. Both the MV methods and the CF ones can narrow the mainlobe compared with the DAS beamformer. However, the speckle performance of the CF is unsatisfactory. The proposed method can obtain a narrow mainlobe and a smooth speckle simultaneously. The deep region performance is also better than the DAS method. Numerical results of the CR and CNR for both cysts are shown in Table 3. The speckle region at the same depth with the cyst is selected for the CR and CNR calculation in order to avoid signal attenuations. From the results, it is shown that the CF can bring an extremely high CR since it can effectively suppress the noise while the CNR suffers a great degradation due to the damaged speckle. On the contrary, MV based methods could bring a higher CNR, but the CR improvement is comparably little. In comparison, the proposed method can significantly increase the CR while the CNR drops a little when compared with the JTR method, but is still at a similar level with the DAS. Sidelobe noises inside the cyst can be effectively removed by the proposed method, which results in a 154% higher CR than the DAS method.

Experimental Phantom Study
The wire target phantom results are shown in Figure 5, which are similar to the simulated results presented before. Both the MV methods and the CF ones can narrow the mainlobe compared with the DAS beamformer. However, the speckle performance of the CF is unsatisfactory. The proposed method can obtain a narrow mainlobe and a smooth speckle simultaneously. The deep region performance is also better than the DAS method.

Experimental Phantom Study
The wire target phantom results are shown in Figure 5, which are similar to the simulated results presented before. Both the MV methods and the CF ones can narrow the mainlobe compared with the DAS beamformer. However, the speckle performance of the CF is unsatisfactory. The proposed method can obtain a narrow mainlobe and a smooth speckle simultaneously. The deep region performance is also better than the DAS method.   Figure 6 gives the lateral variation across the depth of 9 mm. Table 4 provides the statistical results of the FWHM. The proposed JTR-ESBMV GCF obtains a 51% smaller FWHM than the DAS beamformer.  Figure 6 gives the lateral variation across the depth of 9 mm. Table 4 provides the statistical results of the FWHM. The proposed JTR-ESBMV GCF obtains a 51% smaller FWHM than the DAS beamformer.  Results of the anechoic cyst phantom are presented in Figure 7. The advantages of the proposed method mainly occur in better noise reduction, clearer cyst edges, and better-preserved speckle patterns. The CRs and CNRs are given in Table 5 as the evaluation metrics. The experimental results are similar to the simulated ones, in which the JTR-ESBMV GCF achieve the highest CR (111% higher than the DAS one). In addition, the CNR is also a little better than that of the DAS beamformer, but slightly worse than the JTR methods. This indicates that the most prominent advantage of the proposed method is the high contrast. Results of the anechoic cyst phantom are presented in Figure 7. The advantages of the proposed method mainly occur in better noise reduction, clearer cyst edges, and better-preserved speckle patterns. The CRs and CNRs are given in Table 5 as the evaluation metrics. The experimental results are similar to the simulated ones, in which the JTR-ESBMV GCF achieve the highest CR (111% higher than the DAS one). In addition, the CNR is also a little better than that of the DAS beamformer, but slightly worse than the JTR methods. This indicates that the most prominent advantage of the proposed method is the high contrast.    Study   Figures 8 and 9 show in vivo results from two human carotid arteries. The anatomic structure of a blood vessel is similar to an anechoic cyst inside the speckle region. Therefore, the most interesting issues are noises inside the blood vessel and edges of the artery wall. From the results, it could be observed that noises inside the vessel are obvious with the DAS beamformer. Adaptive beamformers can help reduce the noise level while the MV can further preserve the outside-vessel region. The CF result suffers from a low intensity due to the over-suppression of the desired signal. Using the proposed method, the hyperechoic structures can be better distinguished and the noise reduction in anechoic regions is significantly enhanced. Over-suppression is also avoided.  Study   Figures 8 and 9 show in vivo results from two human carotid arteries. The anatomic structure of a blood vessel is similar to an anechoic cyst inside the speckle region. Therefore, the most interesting issues are noises inside the blood vessel and edges of the artery wall. From the results, it could be observed that noises inside the vessel are obvious with the DAS beamformer. Adaptive beamformers can help reduce the noise level while the MV can further preserve the outside-vessel region. The CF result suffers from a low intensity due to the over-suppression of the desired signal. Using the proposed method, the hyperechoic structures can be better distinguished and the noise reduction in anechoic regions is significantly enhanced. Over-suppression is also avoided.

In Vivo
For quantitative assessment, the CRs and CNRs are calculated using the anechoic region inside the blood vessel and the background tissues in Figure 8. Results of different beamformers are given in Table 6. The proposed JTR-ESBMV GCF achieves a high performance in the CR. Statistical results further validate the effectiveness of the proposed method in enhancing the imaging quality.   For quantitative assessment, the CRs and CNRs are calculated using the anechoic region inside the blood vessel and the background tissues in Figure 8. Results of different beamformers are given in Table 6. The proposed JTR-ESBMV GCF achieves a high performance in the CR. Statistical results further validate the effectiveness of the proposed method in enhancing the imaging quality.  For quantitative assessment, the CRs and CNRs are calculated using the anechoic region inside the blood vessel and the background tissues in Figure 8. Results of different beamformers are given in Table 6. The proposed JTR-ESBMV GCF achieves a high performance in the CR. Statistical results further validate the effectiveness of the proposed method in enhancing the imaging quality. The human thyroid results are given in Figure 10. Compared with the DAS result in Figure 10a, both the JTR method and the proposed method can preserve the speckle patterns inside the thyroid, especially at the center of these images. However, the CF result suffers a low intensity and dark artifacts. It is due to the over-suppression of the desired signals caused by a small CF. In comparison with the JTR method, the contrast around the hyperechoic region is enhanced in Figure 10f. The over-suppression could also be effectively avoided. As a result, a better visualization of the thyroid structure can be obtained with the proposed method. The human thyroid results are given in Figure 10. Compared with the DAS result in Figure 10a, both the JTR method and the proposed method can preserve the speckle patterns inside the thyroid, especially at the center of these images. However, the CF result suffers a low intensity and dark artifacts. It is due to the over-suppression of the desired signals caused by a small CF. In comparison with the JTR method, the contrast around the hyperechoic region is enhanced in Figure 10f. The oversuppression could also be effectively avoided. As a result, a better visualization of the thyroid structure can be obtained with the proposed method.  Figure 11 presents in vivo results of a human parotid. Similar to the thyroid results, the hyperechoic structure is more distinguishable with the proposed JTR-ESBMV GCF beamformer. In Figure 11b, the whole image is hard to recognize because of the low intensity. The speckle patterns are severely damaged by the CF. In Figure 11e,f, speckles are well preserved, and the boundary between the hyperechoic region and the speckle is also better defined.  Figure 11 presents in vivo results of a human parotid. Similar to the thyroid results, the hyperechoic structure is more distinguishable with the proposed JTR-ESBMV GCF beamformer. In Figure 11b, the whole image is hard to recognize because of the low intensity. The speckle patterns are severely damaged by the CF. In Figure 11e,f, speckles are well preserved, and the boundary between the hyperechoic region and the speckle is also better defined. All in vivo results demonstrate that the proposed method is effective in enhancing the imaging details and preserving the imaging structure. A better visualization of anatomical structures can be obtained, which could contribute to future medical diagnosis.

Robustness to Sound Velocity Inhomogenetities
An additional experiment is specially conducted to study the influence of sound velocity All in vivo results demonstrate that the proposed method is effective in enhancing the imaging details and preserving the imaging structure. A better visualization of anatomical structures can be obtained, which could contribute to future medical diagnosis.

Robustness to Sound Velocity Inhomogenetities
An additional experiment is specially conducted to study the influence of sound velocity inhomogenetities on different beamformers. When ultrasound waves propagate through different kinds of tissues, such as skins and blood, the sound velocity difference of these tissues could bring calculation errors for the propagation time. In particular, for the PWC, sound velocity errors also result in an inaccurate steering vector. The high performance of conventional beamformers could be degraded as a result [37,38]. Figure 12 shows the simulated PSFs with an 0%, 2%, 5%, and 10% overestimation of the sound velocity using the DAS, the CF, the JTR-MV, the JTR-MV GCF, the JTR-ESBMV, and the JTR-ESBMV GCF. As shown in Figure 12a, the sound velocity errors cause a wide mainlobe and severe sidelobe artifacts. The CF method could effectively suppress the sidelobes, which validates the noise suppressing advantage of the CF. Using the JTR method, sidelobes are partly removed. With the proposed JTR-ESBMV GCF, the sound velocity errors show little influence on the PSF, which can be proved from the comparison between different overestimations in Figure 12f. To conclude, the proposed beamformer achieves a high robustness to sound velocity inhomogeneities.

Discussion
The proposed method is specially designed for the PWC. The MV beamformer is adopted on both apertures to improve the resolution. After obtaining the MV beamformed output matrix, the modified 2-D GCF is used to constrain the final output for the better noise reduction. From the simulations and experiments, the proposed method shows its potential in achieving the high resolution and contrast simultaneously. The innovations of this method can be summarized into three aspects. First, the MV beamforming process is conducted on both apertures. Compared with conventional MV implementations on the receiving aperture, our method fully utilizes the 2-D data information. As a result, the high resolution can be expected. Second, the GCF is specially modified for the 2-D echo data matrix. It is proposed to use the 2-D FFT for the GCF calculation. The lowfrequency component is defined with more accuracy than the original GCF. Thus, the improved GCF could bring better noise reduction and avoid over-suppression of desired signals. Third, the subarray division in the spatial smoothing technique is upgraded into the sub-matrix division, which is applied into both the MV process and the GCF one. Specifically, using the spatial smoothing, the MV helps

Discussion
The proposed method is specially designed for the PWC. The MV beamformer is adopted on both apertures to improve the resolution. After obtaining the MV beamformed output matrix, the modified 2-D GCF is used to constrain the final output for the better noise reduction. From the simulations and experiments, the proposed method shows its potential in achieving the high resolution and contrast simultaneously. The innovations of this method can be summarized into three aspects. First, the MV beamforming process is conducted on both apertures. Compared with conventional MV implementations on the receiving aperture, our method fully utilizes the 2-D data information. As a result, the high resolution can be expected. Second, the GCF is specially modified for the 2-D echo data matrix. It is proposed to use the 2-D FFT for the GCF calculation. The low-frequency component is defined with more accuracy than the original GCF. Thus, the improved GCF could bring better noise reduction and avoid over-suppression of desired signals. Third, the subarray division in the spatial smoothing technique is upgraded into the sub-matrix division, which is applied into both the MV process and the GCF one. Specifically, using the spatial smoothing, the MV helps obtain the beamformed outputs of different sub-matrices. After that, the 2-D GCF is computed using the MV output rather than the original echo data matrix. Since the sub-matrices are overlapped, the robustness is therefore enhanced and a smooth speckle performance can be obtained. In consideration of these factors, the high performance of the proposed method is reasonably convincing.
The proposed method has a comparably lower CC than conventional adaptive beamformers. For each PWC with N emissions, it only needs to calculate two MV weights and one GCF while single frame methods require N MV weights and N GCFs [16,24]. Due to the application of the FFT, the GCF calculation has a much lesser computational amount than the MV process and can be neglected. Therefore, the computational load mainly comes from the dual aperture MV process. For receiving aperture beamforming, the computation for ESBMV weights has a complexity of O L 3 1 while for transmitting aperture, the number becomes O L 3 2 [19,21,36]. Specifically, for the L 1 × L 1 size sub-matrix in the receiving aperture, the matrix inversion and eigen-decomposition require 2/3L 3 1 and 21L 3 1 floating operations, respectively. Likewise, the same calculation can be done for the transmitting aperture. For the GCF calculation, in consideration of the FFT technology, the 2-D Fourier transform of the (N − L 2 + 1) × (M − L 1 + 1) size matrix requires (N − L 2 + 1) log 2 (N − L 2 + 1) + (M − L 1 + 1) log 2 (M − L 1 + 1) operations. The total computational amount is 2/3L 3 1 + 21L 3 1 + 2/3L 3 2 + 21L 3 2 + (N − L 2 + 1) log 2 (N − L 2 + 1) + (M − L 1 + 1) log 2 (M − L 1 + 1), only slightly higher than that of the JTR method. In addition, it can be inferred that there is a positive correlation between the calculation amount and the subarray size. It accords with the statement presented earlier that a larger subarray size could bring better noise reduction and lower sidelobes at a cost of more calculations. In consideration of these factors, the proposed method is promising for ultrafast ultrasound imaging.
Thanks to the wide application of the graphics processing unit (GPU), the real-time MV implementation could be realized [39,40]. However, the static imaging and dynamic videos are quite different. Specifically, in practical situations, tissue motions and environment disturbance could severely degrade the imaging quality. In addition, whether the proposed method introduces frame-to-frame artefacts is worth investigation. Therefore, it would be valuable to study the video data of a real-time implementation. A comparison between the dynamic and static imaging is also necessary. Future studies will focus on the real-time in vivo application of the proposed method. Related techniques in audio processing will also be taken into consideration [41].

Conclusions
In this paper, we proposed a novel beamformer for PWC, which integrates the 2-D MV method with the GCF. Both the MV and the GCF were modified for the PWC imaging mode to maintain the high resolution and contrast. Results of simulations, phantom experiments, and in vivo studies validate that the proposed method was effective in enhancing the imaging quality. A high contrast and a high robustness against the sound velocity error could be obtained. In conclusion, the proposed method could be promising for maintaining the satisfactory imaging quality for ultrafast medical ultrasound imaging. Future works will focus on the real-time application of the proposed method, especially in in vivo situations. The tissue motion compensation technique will be taken into consideration to improve the dynamic imaging quality.

Conflicts of Interest:
The authors declare no conflict of interest.