A 2D phase zero algorithm for estimation of displacement in ultrasound elastography based on the optimization of initial value for iteration

ABSTRACT Ultrasound elastography (UE) is an imaging technique based on tissues’ elasticity. UE based on gradient typically assumes that movement only occurs along the propagation direction of an ultrasound beam. However, when the target tissue is compressed by a transducer, lateral displacement of the scattering body occurs due to uniform hardness within the tissue. Before and after compression, the signal cross-correlations within a specific window along the scan lines may be very low. This is highly likely to result in large errors or even failed calculation of cross-correlated phase angles within this window. We proposed a modified 2D phase zero algorithm for the estimation of displacement. In this method, the initial value with threshold for iteration was optimized in order to decrease the errors result from lateral motion. Then, to verify the effectiveness of this algorithm, a contrast experiment on this new and conventional method was designed. The experiment indicated that this algorithm could effectively prevent failed calculation of elasticity due to lateral displacement of the scattering body. It could stably and efficiently reduce the errors in displacement estimation and enable a more accurate imaging of the target hard mass, like signal to noise of elastography () and contrast to noise of elastography ().


Introduction
Ultrasound elastography (UE) is an emerging medical imaging modality (E-mode) which provides a graphical representation of tissues' elasticity, either directly or indirectly [1]. Thus physicians can make pathological and physiological assessments of the target tissues based on the correlation between pathological changes and hardness/elasticity changes of the tissues. UE is particularly useful for early detection and diagnosis of the lesions.
UE usually consists of three steps. These steps are the following: displacement estimation, strain estimation and strain display. Displacement can be estimated based on the maximum of the cross-correlation function in the time domain or by finding the phase root in the crosscorrelation between the pre-and post-compression complex base-band echoes. Alternatively, instead of using the cross-correlation function, displacement can also be estimated by finding the absolute value or sum of squares of the phase difference [2][3][4][5].
Algorithms relating to the phase domain can produce an estimation of displacement with high speed and low cost. The most common phase-domain algorithm is the phase zero algorithm proposed by Ophir's team [2]. This algorithm assumes that the sampling interval between the two frames is very small before and after compression. It also assumes that the compressed hard tissue makes a rigid parallel motion. In other words, only axial displacement occurs in the two frames. Then the displacement is estimated through the cross-correlation of signals of the two frames within the 1D vertical window. However, the hard mass is not an ideally rigid body. It also has heterogeneous hardness and makes horizontal deformation and motion. Therefore, displacement not only occurs in the window of the two scan lines before and after compression, but also occurs in the window around the two scan lines. Consequently, the estimated displacement containing error will be used as the initial value of iteration for the next window.
Addressing the cross-correlation problem between the 1D windows, Jiang et al. [6] proposed a 2D phase zero algorithm for displacement estimation. Integrating the information of several scan lines in this algorithm, it extended the 1D windows into 2D windows. This enhanced the cross-correlations between the windows and reduced the influence from side slipping of the 1D axial signals. By this means, the accuracy and robustness of the algorithm were greatly improved. However, this 2D extension only enhanced the cross-correlations between the windows. Under large compression, tissue makes lateral displacement as a mass, which leads to calculation errors even with the 2D phase zero algorithm. Given the high dependence on the initial value for iteration, any error in the initial value using the phase zero algorithm will be magnified in the following windows. This finally results in failed estimation of the axial displacement. Correction of the initial value is one way to prevent error magnification in the phase zero algorithm.
Based on the 2D phase zero algorithm, Zhang et al. [7] proposed a method of initial value correction based on statistics in the horizontal windows (LDS-PZA), which utilized the high cross-correlations in axial motion along adjacent scan lines. Axial displacement between the scan lines was considered very small. Furthermore, the calculated results in several adjacent axial windows were analyzed statistically to remove the wrong displacement values. The remaining values were averaged to get the initial value for iteration in the following axial window. However, this method does not judge whether rigid deformation exists or not. Neither does it make a careful differentiation of all calculated initial values in the windows. By directly removing the smallest value, this algorithm fails to eliminate other wrong values close to the smallest value and mistakenly eliminates the valid values.
To ensure a good estimation of displacement, many researchers made a lot of optimization methods. In 2015, Bal and Imperiale [8] proposed a reconstructed procedure for the spatial variation of the elastic displacement vector structure in biological tissue. An optimization method has been proposed based on the monogenic features combined with robust filtering technique by Slimi et al. [9]. Zhang et al. [10] proposed a 2D ultrasonic elastography based on beam steering and iterative correction to improve the quality of imaging.
With the rapid development of computer technologies [11], we proposed a 2D phase zero algorithm for displacement estimation based on the optimization of the initial value for iteration. Our algorithm is innovative in that it includes a reliable differentiation of the initial values by setting the threshold.

Materials and methods
1d and 2d phase zero algorithm Let the signal of tissue compression be x 1 .The echo signal x 2 is considered a time-shifted signal of x 1 . Supposing that within a certain range centering around time t, the signal only undergoes rigid parallel motion but no stretching, x 2 is written as.
where t is the displacement (time shifting) to be calculated. Cross-correlation between x 1 and x 2 is represented by h x 1 ; x 2 i ðtÞ, which is itself equivalent to an autocorrelation of x 1 : Since this autocorrelation has the maximum value at t ¼ 0, the cross-correlation has the maximum at t ¼ Àt. Cross-correlated phase ' t ð Þ also exists at t ¼ Àt, for x 1þ t ð Þ and x 2þ t ð Þ, the analytical signals of x 1 and x 2 , respectively.
By finding the solution to ' t ð Þ ¼ 0, the displacement t can be estimated. ' t ð Þ is approximated as a linear representation using the central frequency of the transducer v 0 as the slope. The solution to ' t ð Þ ¼ 0 can be found by Newton's iteration method: However, the actual sampled signals are discrete data. Linear interpolation on the baseband In-phase and Quadrature (IQ) data is performed. Demodulation is the process of removing the carrier frequency and reconstructing the fundamental signals, expressed as follows: where v m is the demodulated frequency, usually replaced by the central frequency of transducer v 0 . Therefore, the cross-correlation of the analytical signals is expressed as follows: From Equations 4-6, the algorithm for time delay estimation in the discrete domain is expressed as follows: where t k;l is the time shifting in the k-th window after l iterations; T w is the window size used for cross-correlation calculation; Ã is the complex number conjugate.
Due to non-uniform hardness of the tissue in the 1D phase zero algorithm, the scattering body will undergo lateral displacement. As a result, before and after compression, the signal cross-correlation within a specific window along the scan line is very low. This will lead to calculation errors of the cross-correlated phase within this window or even total failure. Zhang et al. [7] propose the 2D phase zero algorithm: Ás 2b x; t À t x;yÀ1 2 dt where x is the horizontal direction and y the vertical direction; t x;y is the estimated time delay of the x-th scan line at the y-th position; T w is the size of the vertical window; n is the size of the newly added horizontal window.

Optimization of initial value for iteration based on statistics in horizontal windows
The 2D phase zero algorithm involves the calculation of the horizontal windows to reduce the influence from the horizontal deformation or displacement of the scattering body. It is because of this additional step that there is a wrong estimate of the initial value for other reasons. Zhang et al. [7] proposed a statistical method for initial values to remove the wrong values. The 2D phase zero algorithm can be represented with Equation (9): t x;yÀ1 ¼ avg t xÀ T l 2 ;yÀ1 :::t xþ T l 2 ;yÀ1 À min abs t xÀ T l 2 ;yÀ1 :::t xþ T l where min denotes the finding of the minimal value; T 1 and T a represent the horizontal and vertical window, respectively.

2D phase zero algorithm based on optimization of the initial value for iteration
However, the external excitation applied to the tissue was attenuated over a larger depth. Therefore, the horizontal deformation and displacement of tissue become less likely at a larger depth, which corresponds to reduced probability of a wrong displacement estimate. Since the displacement between the scan lines increases with depth, the calculation result will be distorted considerably if the horizontal windows used for statistics in the far field is too large. The correction method described in Equation 9 does not provide an effective discrimination between the initial values in the statistical windows. Instead, it directly removes the minimal value. This results in two situations: one is the failure to remove other wrong values close to the minimal value; the other is the removal of valid values. In order to solve this problem, we proposed a 2D phase zero algorithm based on optimization of the initial value for iteration.
Equation (9) is modified to: t x;yÀ1 ¼ avg t xÀ T l 2 ;yÀ1 :::t xþ T l 2 ;yÀ1 À f re abs t xÀ T l 2 ;yÀ1 :::t xþ T l where a is a threshold above 0 but below 1; f re denotes the solving of the relative error, i.e. percentage of mean offset of statistics, found by the following: f re ðmÞ ¼ avg t 1 ; t 2 :::t n ð ÞÀabs t m ð Þ avg t 1 ; t 2 :::t n ð Þ The procedures of the 2D phase zero algorithm based on the optimization of the initial value for iteration are presented as follows ( Figure 1): 1) In-phase and Quadrature (IQ) data were divided into M £ N windows based on the preset parameters of 2D windows (window size and overlap rate). The initial value for iteration in the first row of windows was set to 0; 2) Using the 2D phase zero algorithm, the displacement of the m-th (m = 1,2,…M) axial window was calculated and recorded for each column; 3) Slide the statistical windows. The statistics of the m-th axial window were also obtained. It was then determined whether the relative error was above the threshold a. If yes, remove the values. Also, the average was taken as the displacement for the corresponding window; 4) The statistics of windows in all columns were taken as the initial values of the m + 1-th axial window; 5) Step 2 was repeated until the calculation was performed for all axial windows.
The size of the horizontal windows changes with depth. That is, a larger window is used for the near field, and a smaller window is used for the far field, so as to ensure the accuracy of displacement estimate. The size of horizontal window T 1 was calculated by Equation (12): where S is the number of sampling points along the scan lines; s th is the vertical position where the displacement is estimated; T L 0 is the initial size of the horizontal window.

Simulation platform
The phantom in the experiment was designed and used as a test object in [12]. The uniform phantom consists of 1.5 £ 10 5 scatterers (fully developed) distributed uniformly throughout the dimensions of 40 mm £ 30 mm £ 5 mm with scattering strengths Gaussianly distributed. The IQ signals with a 128-channel transducer and 7.5 MHZ central frequency of each frame were simulated using Field II, which is a computer program designed by Jensen JA and described in [13]. In this simulation experiment, the threshold value was set to 7.3%.

Results and analysis
Two algorithms, LDS-PZA and MLDS-PZA, were applied respectively to estimate displacement under different compression ratios in elastography. Figure 2 is the comparison of the axial strain imaging using the two algorithms, respectively. The images obtained using the two algorithms under different compression ratios were analyzed: In Figures 2(a-b) and 3(a-b), when the compression ratio was 2% to 3%, both of the two algorithms achieved satisfactory results, with an effective control of errors in the estimated displacement; In Figures 2(c) and 3(c), when the compression ratio was 4%, errors appeared at the margins of the strain diagrams for both of the two algorithms, which were small in area and confined to the two sides; In Figures 2(d) and 3(d), when the compression ratio was 4.5%, the LDS-PZA algorithm led to an abrupt increase of errors. This is because of several displacements with smaller absolute values in the 2D window. However, with the LDS-PZA algorithm, only one minimal value was removed, leading to significant errors and an accumulation of errors. As a result, errors not only appeared at the margins, but they also affected the entire diagram, indicating total failure. In contrast, the MLDS-PZA algorithm achieved a much more accurate estimate. Continuous errors only occurred locally on the two sides, while the displacements were clearly seen for the target region. In Figures 2(e) and 3(e), when the compression ratio was 5%, the imaging result was much improved using the LDS-PZA algorithm as compared with the result at 4.5% compression ratio. However, the dark vertical strips still covered the target region. The result with the MLDS-PZA algorithm was more satisfactory, especially in the target region.
Finally, in order to quantitatively analyze the image quality of the strain image, the SNR e and CNR e were calculated based on LDS-PZA and MLDS-PZA, as shown in Figure 4.
In Figure 4(a), the result images based on those two algorithms were able to maintain a high SNR e in the case of compression ratios less than 3%. But when the com- pression ratio is in the range of 4%-5%, those results were affected by the increase of edge error results, and the value of SNR e decreases rapidly. However, MLDZ-PZA is slower than LDS-PZA in the descent rate of SNR e .
In Figure 4(b), the result images based on those two algorithms were able to maintain a high CNR e in the case of compression ratios less than 4%, which means that the target block is well protected, and MLDS-PZA has a smaller CNR e than LDS-PZA. When the compression ratio is in the range of 4.5%-5%, the CNR e of the two algotithms decreases rapidly, but MLDS-PZA is more stable.
Combined with the above results, MLDS-PZA uses a threshold to distinguish the initial values, which can remove other error data that are the same or close to the minimum, also retaining valid data. This makes MLDS-PZA more stable and more effective in reducing the displacement estimation errors in high compression ratios, and more accurately imaging the target hard block region.
However, the MLDS-PZA algorithm is limited in the choice of threshold. Since the threshold value is related to the compression ratio, the compression ratio should be one of the parameters considered for setting the threshold. The compression ratio is usually unknown in practice. An additional model may be needed to determine the compression ratio. For example, the block matching algorithm could be used. Therefore, estimation of the compression ratio will be addressed in future investigation.

Conclusions
To improve the quality of elastography, we analyzed the defects in the conventional 2D phase zero algorithm involving the statistics in horizontal windows. Finally, we proposed a modified algorithm based on the optimization of the initial value for iteration, MLDS-PZA. Our algorithm is innovative in that it includes a reliable differentiation of the initial values by setting a threshold.
This makes MLDS-PZA more stable and more effective in reducing the displacement estimation errors. The results of the simulation experiment indicated that, when compared with LDS-PZA, the MLDS-PZA algorithm achieved better differentiation between the initial values. Also, through an effective correction of the estimated displacements, it had higher stability. Even at high compression ratios, the MLDS-PZA algorithm was still capable of producing the strain diagram for the target region. Although further improvements are needed, the algorithm presented here is promising in that it has greatly enhanced the accuracy of imaging and is helpful for clinical diagnosis and treatment.