1 Introduction

In order to improve the level of medical treatment, reduce medical costs and meet the basic medical needs of the masses, telemedicine has been launched as a new medical service. At present, telemedicine technology has realized the real-time communication of voice and high-definition image using high-speed network, which provides a broader development space for modern medicine [1,2,3]. And national health authorities are focusing on E-health services such as E-health cards, health portals, health information networks, electronic medical records, telemedicine services, and personal wear and portable communication systems for patient monitoring and support [4,5,6]. Among them, independently managed biomedical/medical stakeholders (such as hospitals, suppliers, and patients) will adopt blockchain technology for applications that want to collaborate without ceded control to centrally managed intermediaries [7,8,9,10,11].

Today, COVID-19 leads to reduced lung function in patients, leading to death [12]. As an indispensable part of telemedicine and pulmonary ultrasonography, real time ultrasound elastography has attracted much attention. It expands the connotation of ultrasonic diagnosis theory and the scope of ultrasonic diagnosis, makes up for the deficiency of traditional ultrasonic imaging, can show the lesions more vividly, differentiate the nature of the lesions [13, 14]. It reflects the physical information of soft tissue through the image, discriminates the elastic size of pathological tissue, and thus deduces the possibility of certain pathological changes. Because the traditional ultrasound imaging method provides only a reflection of the acoustic impedance of soft tissue [15, 16]. Traditional ultrasound imaging methods have limited recognition ability, especially for soft tissue with small acoustic impedance difference. Moreover, the detection of deep body tissues is more likely to be interfered by various factors, such as shielding, light or shadow, and the patient’s physical condition. Ultrasound elastography technology can address this deficiency. Ultrasound elastography has become one of the research topic of interest in the ultrasonic imagery field both at home and abroad [17,18,19,20,21], and various improved ultrasonic imaging methods have emerged, and aroused much research interest.

For example, Ophir J.’s method [22] can create a slight deformation of the tissue by applying a little pressure and collect an ultrasonic signal through an ultrasonic device. Using a time-domain cross-correlation to compare the two sets of signals, we can obtain the displacement of each point in the tissue and calculate the displacement difference to obtain the tissue strain. The reciprocal of the strain is used to appropriately represent the physical attributes of the tissue. The precision of cross-correlation in time domain is mainly affected by ultrasonic frequency, sampling rate and signal-to-noise ratio. Therefore, when the SNR is high, it is more sensitive to the sampling rate [23, 24]. In O’Donnell M.’s method, the signal envelope cross-correlation method is used to calculate the tissue displacement [25]. In this method, a Hilbert transform is performed on the collected RF to obtain the analytical signal, and the analytical signal is cross-correlated. we can obtaine the displacement data of the tissue by analyzing the phase of the analytical signal cross-correlation function. Due to the periodic characteristics of the phase, this method can only be used when the maximum displacement is less than a 1/4 wavelength. When the displacement is greater than a 1/4 wavelength, phase cancellation will occur and produce an incorrect displacement and method failure.

To address this problem, a CAM (Combined Autocorrelation Method) [26, 27] is proposed, which expanded the application scope of the method. But the CAM requires a large amount of computation and is greatly affected by noise. The block-matching method is widely used in video compression and video tracking applications [28]. It divides the target frame into blocks of the same size. The best matching position is found in the reference frame (usually the adjacent frame of the target frame) for each block. There are many measures of matching degree, including SAD (Sum of Absolute Differences) [29], MSD (Mean Square Differences) [30] and NCC (Normalized Cross Correlation) [31]. The best matching location search methods include FS (Full Search) [32], TSS (Three-Step Search), NTSS (New Three-Step Search) [33], and 4SS (Four-Step Search) [34]. To narrow the search scope, the search can be carried out around the best matching position that has been calculated for the adjacent blocks. However, because this search method based on prior knowledge, it can lead to error accumulation and obtain intolerable error results. In this paper, the block matching method and zero phase method are improved and combined. We propose a novel displacement estimation method for real time tissue ultrasound elastography and achieve good experimental results. The former provides lead information so that the latter can get good results even in the case of large displacement. That is, a quality-guided block matching method is adopted to improve the computing speed and avoid the accumulation of errors.

Therefore, the organizational structure of this paper are: Section 1 introduces the research background of ultrasound elastography in the mobile telemedicine. Section 2 presents the quality-guided block matching method and the phase-zero search method and introduces our method framework. Section 3 introduces two comparative methods including CAM and time delay estimation (TDE). We make some comparitive experiments in Section 4. Section 5 concludes this paper and outlook the future work.

2 Related work

2.1 Quality-guided block matching

In the field of digital image processing, the sum of absolute differences (SAD) is often used to measure the similarity between image blocks. It first gets every pixel in the original block. Then calculate the absolute difference with the corresponding pixel in the block for comparison. Finally, these differences are added together to create a simple block similarity measure [35, 36]. At present, the sum of absolute differences is often used for object recognition, motion estimation for video compression, generation of disparity maps for stereo images, and so on [37]. Because of its simple calculation, convenient operation and fast calculation speed, and considering the simplest possible metric of each pixel in the block, it is very effective for a wide motion search of many different blocks. In addition, SAD can analyze each pixel individually and is easy to parallelize processing, so it is easy to implement.

In our method, SAD as a measurement of matching degree is defined as :

$$ SAD = \sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}|X(i,j)-Y(i,j)| $$
(1)

where, m × n is the size of the block, X is the block to be computed in the target frame and Y is the block of the reference frame.

Corresponding to the reference frames, the best matchable position in block X(k,s) is:

$$ \begin{array}{ll} Y(k+p^{\prime},s+q^{\prime}) = \arg\min(\sum\limits_{i=k}^{k+m}\sum\limits_{j=s}^{s+n}|X(i,j) -\\Y(i+p,j+q)||p\in(-x,x),q\in(-y,y)) \end{array} $$
(2)

where, (−x,x) is the transverse searching scope and (−y,y) is the vertical searching scope. Meanwhile, we take the displacement of the points in block X as a vector \((p^{\prime },q^{\prime })\).

On the basis of the traditional block searchin, which computes the displacement of every block by either row order or column order, using the result of the neighboring block of the current block Xs a reference computes the displacement of the current block. Due to the effect of noise, errors can accumulate continuously and even result in a large error if there is a displacement error in one block [38, 39].

Block matching is a common method in image denoising and motion estimation. By matching the query block with the adjacent image block, K blocks with the closest distance from the query block are found out from these adjacent blocks. The quality-guided block matching method [40] does not calculate the blocks in the target frame in the order of row and column. it calculates the displacement of its neighbor block by referring to the block with the highest matching degree at a given point in time. Using this method can ensure that the more matchable block gets calculated first, and the block that has more noise and easily results in error is computed last, thus avoiding the transmission and accumulation of errors.

The steps of quality-guided block matching method are as follows:

  1. step 1.

    Select the starting block X and put it into the set S to calculate the displacement and matching degree of block X;

  2. step 2.

    Find the most matchable block in set S and take it as the current block. If there is more than one candidate, then randomly select from the candidates;

  3. step 3.

    Compute the displacement and matching degree of the current block;

  4. step 4.

    If there are neighboring blocks of the current block that have not been computed, go to Step 5; otherwise, go to Step 7;

  5. step 5.

    If the neighboring block chosen in the last step was not in set S, initialize the neighboring block with the displacement and matching degree of the current block, put it into set S and return to Step 4. If it was in set S, go to Step 6;

  6. step 6.

    If the matching degree of neighboring block is higher than the current block, return to Step 4; otherwise, update the displacement and matching degree of the neighboring block with that of the current block and return to Step 4;

  7. step 7.

    Remove the current block from set S;

  8. step 8.

    If set S is not empty, return to Step 2, otherwise the steps are finished.

2.2 Phase-zero search

The phase-zero search was first proposed by Cabot in 1981. The main principle is to calculate the analytic form of the signal of the two frames before and after press-down, and calculate the position of the phase crossing zero by using the analytic signal cross-correlation to get the maximum position of the cross-correlation function. In other words, the number of relations between each other will get real digital signal structure for complex analytical signal, the analytical signal phase has the following properties: in relationship to each other for digital signal maximum phase near zero, corresponding to the precise relationship between each other several continuous signal of the maximum point, through the precise the size of the maximum point position deviation disposal, can get relatively precise time delay signal [41, 42].

In general, we do not consider the relationship between hardware devices. We can use the software to calculate the zero position, its principle is to echo signals structural analysis, and seeking cross-correlation function of complex signal in time domain to calculate the cross-correlation function. When you have the maximum t = 0, for analytical signal cross-correlation function has a phase of zero. Finally, by Newton’s method to calculate the time delay.

In fact, the phase-zero search first calculates the analytical form of the two frames of signals before and after palpation, and then, It uses the cross-correlation of analytic signal to calculate the position of zero crossing zero, that is, the position of maximum value of cross-correlation function. The bandpass of the two signals x1(t) and x2(t) are described as:

$$ \widetilde{x}_{1} = A(t - \tau_{1})e^{-i\omega_{0}\tau_{1}}, \widetilde{x}_{2} = A(t - \tau_{2})e^{-i\omega_{0}\tau_{2}} $$
(3)

The complex cross-correlation functions of the two signals are described as:

$$ \widetilde{C}(t)= \frac{1}{T}{{\int}_{0}^{T}}\widetilde{x}_{1}(\tau)\widetilde{x}^{*}_{2}(t+\tau)d\tau $$
(4)

When t is 0, we can obtain the simplified form of the above expression:

$$ \widetilde{C}(0) = {\Gamma}_{AA}(\tau)e^{-i\omega_{0}\tau_{1}} $$
(5)

where, ΓAA(τ) is an autocorrelation function of the signal envelope.

The time delay of signals x1(t) and x2(t) can be described with the above phase-zero search:

$$ \tau_{BB} = \frac{\phi(0)}{\omega(0)} = \frac{\tan^{-1}(\frac{Im(\widetilde{C}(0))}{Re(\widetilde{C}(0))})}{\omega(0)} $$
(6)

The displacement is:

$$ d = \tau_{BB} \cdot c $$
(7)

where, c is the propagation speed of ultrasound in the medium.

However, the disadvantage of this algorithm is that the displacement accuracy is proportional to the number of iterations. The higher the displacement accuracy, the more iterations. That is to say, the method is only applicable to the case of small wavelength, can not be used for large displacement.

2.3 Method framework

The method framework for this paper is shown in Fig. 1. The phase-zero search can obtain equal accuracy with oversampling by a lower calculation and cause phase cancellation, resulting in a miscalculation of time delay when the displacement is more than a half phase. Combining block-matching with a phase-zero search can obtain an accurate displacement field after first computing a rough displacement with block-matching that limits the error no more than a half phase.

Fig. 1
figure 1

The algorithm framework

Full size image

To clearly reflect the physical characteristic of the tissue, we compute the strain field from the displacement field with the gradient method. Specifically, we first used RF signal to get a rough displacement estimation, and then analytical signal and rough displacement estimation to get an accurate estimation. Because if accurate displacement is obtained directly, this result may not be accurate. If rough displacement estimation is applied again, it will have a relatively good effect. Secondly, the RF signal before and after pressure are measured, and analytical signal is obtained by Hilbert Transform and Block-Matching. Rough displacement estimation and analytical signal were used to calculate accurate displacement by the phase-zero search method. Finally, the least squares method is used to get distribution of strain.

3 Comparative approach

3.1 Combined autocorrelation method

The CAM [26, 27] and the traditional TDE [43] are classical and popular methods of the current displacement estimation methods. Therefore, this paper uses these methods as a comparative experiment. Many elastography methods based on biological tissues use static tissue compression methods, which actually measure the strain distribution in the body caused by tissue compression or relaxation. Mechanical methods using stepper motors have been used to compress tissue axially. However, it is very simple and desirable to operate the ultrasound probe with bare hands. Generally speaking, we must be robust to the non-axial movement of the probe on the surface when we compress with our hands. Moreover, for stable measurement that does not depend on compression speed and compression, a large dynamic strain range is also required. In order to satisfy these conditions at the same time, CAM is proposed [44].

The CAM produces elasticity images with high-speed processing and high accuracy, and achieves a wide dynamicrange for strain estimation by combining two-step processing.

  1. step 1.

    Search for a rough estimate of the maximum envelope correlation. The key feature of this method is to search only on the grid points of the 1/2 wavelengthinterval in the axial direction and the interval of scan linesin the lateral direction because the purpose is to detect azone without phase aliasing in a large dynamic range. Thus, an improvement in the processing speed is attainedby this method. Here, the 1/2 wavelength interval in the axial direction enables us to optimize computing efficiency without phase aliasing. Moreover, the CAM is robust tosideslip and suited to freehand compression by implement-ing a 2-D search in lateral directions.

  2. step 2.

    The unpacking phase obtained in the first step is used for better estimation. The time shift, the axial displacement, and the lateral displacement, at a measurement point are given respectively. Finally, the axialstrain distribution can be obtained by differentiating theaxial displacement distribution spatially.

To maintain simplicity, we consider only the axial displacement of the tissue (that is, along the direction of the beam). And from the perspective of envelope correlation, the first step is similar to the speckle tracing algorithm. The first and second steps can occur simultaneously through autocorrelation processing, and the envelope correlation coefficient is computed only a few times. So the processing speed is faster than the speckle tracking, and the correlation coefficient is calculated by changing the time shift many times. In practice, as shown in Fig. 2, by using multiple autocorrelation processing units in parallel, it should be possible to perform real-time processing similar to traditional Doppler methods.

Fig. 2
figure 2

Combined autocorrelation method

Full size image

3.2 Time delay estimation

To calculate the time delay of a reference signal and a contrast signal in a period of time, time delay estimation is proposed. TDE has been widely used because of its low complexity, simple operation, good real-time performance and high acoustic positioning accuracy. As shown in Fig. 3, which shows the time shift exists between the reference signal and the contrast signal. The present TDE method has good performance under high Signal to Noise Ratio (SNR). However, with the decrease of signal correlation, the performance decreases, so the robustness is poor [45, 46].

Fig. 3
figure 3

Reference signal and delay signal

Full size image

Because the tissue vibration is small and its orders of magnitude are usually measured in microns, the sampling frequency of common ultrasonic equipment cannot meet the requirement. Therefore, a signal interpolation algorithm must be used to improve the estimation accuracy of displacement. The interpolation algorithm is an significance function approximation method [47]. The spline interpolation finds a set of fitting polynomials according to the existing data points. In the fitting process, it is a common method to use polynomial to fit the curve of adjacent data points in signal processing.

In signal processing, the mutual transformation of continuous and discrete signals is a basic task, and the spline interpolation is the most suitable method. Schoenberg proposed the theoretical basis of spline interpolation and introduced the B-spline curve. Spline curves are represented by piecewise polynomials and connected together smoothly. Connection points are called nodes. For an n-order spline interpolation, each segment of its polynomial is an n-order, and each segment of the curve requires n + 1 coefficients. The curvature of a cubic spline curve is the smallest, so it is mostly used in practice.

The cubic spline interpolation is typically used in signal processing because of the balance between computation and accuracy. With spline interpolation, the discrete signal can be expressed in a continuous form as a polynomial. Then, the original signal can be interpolated by increasing the sampling frequency or by directly applying the polynomial coefficient to the delay calculation [48].

4 Experiments

A two-dimensional soft tissue model was created using Abaqus software, and this model contained circular tissue that was four times harder than the surrounding tissue. Then, a 2% deformation is generated in the model, and a simulated RF signal is generated by the consistent displacement generated by the deformation. In this set of simulation data, the central frequency of the signal is 5 MHz, the sampling rate is 20 MHz, and the maximum displacement is greater than twice the length of the ultrasonic wave. Our displacement estimation method is compared with the CAM and the traditional TDE. These two methods are often applied to real time tissue ultrasound elastography, but the accuracy and time efficiency of the algorithm still needs to be improved. Therefore, we conducted two experiments: estimation of longitudinal displacement and strain estimation.

4.1 Estimation of longitudinal displacement

Firstly, we estimate the longitudinal displacement. In the TDE method, RF signals are initially oversampled. When looking for the peak of the cross-correlation function, the parabolic interpolation method is adopted. Figure 4a, b and c show the displacement field calculated by the TDE method, the CAM method and our method, respectively, and Fig. 4d shows the comparison between the displacement distribution of the three methods on a longitudinal line and the real displacement distribution.

Fig. 4
figure 4

Displacement calculation comparison of three methods

Full size image

Obviously, the approach we proposed is closer to the real displacement distribution because the method we proposed takes the displacement of the neighboring points as the prior information when calculating the displacement of a certain point, which not only improves the accuracy but also greatly reduces the search scope, thus reducing the amount of calculation and speeding up the calculation speed. To obtain the same result as the method in this paper, the TDE needs to carry out oversampling with a very large amount of computation. There are also many redundant cross-correlation calculations in the CAM, and the half wavelength limitation cannot be broken.

To quantitatively analyze the accuracy of the three methods, the displacement difference of the calculation results is defined as:

$$ \sigma = \sqrt{\frac{{\sum}_{i=1}^{n}(result_{i} - true_{i})^{2}}{n}} $$
(8)

where, result is the displacement field calculated, and true is the error of the three methods (TDE, CAM, our method), n is the number of samples.

After calculating the displacement field error in three methods, the results were as shown in Table 1 (0.1156, 0.0798, 0.0657 respectively). Obviously, the displacement field error of this method is the smallest, which indicates that the calculation result of this method is more accurate. Therefore, the method in this paper had the lowest displacement field error and higher accuracy. Figure 4 also shows that our proposed method is closer to the actual displacement distribution and has a better effect.

Table 1 Displacement field error
Full size table

4.2 Strain estimation

In addition, to clearly indicate the physical characteristic of the tissue, we compute the tissue strain on the basis of the gained displacement field with the least square method, namely, the tissue strain is represented by a gradient. Figure 5a, b, c separately represent the results of the TDE, the CAM and the method we proposed. As seen from the figure, the result obtained by the method in this paper is smoother, and the diseased tissue is easier to distinguish. Figure 5d shows the strain distribution at the midline position.

Fig. 5
figure 5

Strain comparison of three methods

Full size image

A comparative transport rate (CNR) is adopted to quantitatively analyze the lesion degree resolution of the tissue of each method, as follows:

$$ CNR = \frac{m_{b}-m_{t}}{\sqrt{({\sigma_{b}^{2}}+{\sigma_{t}^{2}})/2}} $$
(9)

where, mb,σb are the mean and standard deviation of the background, and mt,σt are the mean and standard deviation of the target. The CNR of there methods (TDE, CAM, our method) are shown in Table 2.

Table 2 A comparative transport rate
Full size table

Through calculating the comparative transport rates of the three methods, the results are shown in Table 1, which are respectively 7.67, 26.08 and 157.06. In general, the greater the comparative transport rate of the method, the faster and more effective the calculation results of the method are. Therefore, the comparative transport rate of the method we proposed is the highest and the method is more efficient. Figure 5 also shows that our proposed method is more effective and stable.

5 Conclusion and future work

In this paper, we propose a novel displacement estimation method for real time tissue ultrasound elastography that combines the block-matching method and the phase-zero search method. The experimental results show that the method we proposed can accurately and efficiently calculate the displacement field and effectively solve the error transmission problem when using prior information.

In the future, we will continue to optimize and improve the displacement estimation method. For example, in the quality-guided block matching method, the determination of block size can be adaptive. In the phase-zero search method, the results can be calculated by an iterative method to make the results more accurate. Moreover, with the development of computer software and hardware, it can be used to accelerate the method to meet real-time requirements. In addition, with the widespread application of mobile telemedicine, large hospitals or specialized medical centers will take advantage of medical technology and medical equipment to diagnose, treat and consult patients in remote areas.