1. Introduction
According to the World Health Organization research data, it shows that one of the three major diseases that threaten human life and health is cardiovascular and cerebrovascular disease [1,2]. In the diagnosis of cardiovascular disease, blood volume flow is a very important health indicator [3]. It is natural that the measurement of volume flow is beneficial to the heart, carotid, renal and visceral circulation, as well as to the aorta and cardiac circulation [4-6]. It is an important method to diagnose the various diseases with the ultrasonic Doppler spectrums to detect the blood volume flow of the related vessels [7,8]. Therefore, the ability to accurately and quickly monitor blood volume flow is very essential for clinicians [9], which in clinical diagnosis is of great significance [10].
The general method for estimating blood volume flow using ultrasound is calculating the integral of the measured spatial mean velocity within the vessel cross-sectional area [11-13]. Boote [14] proposed a method estimating volume flow based on a single gate, he estimated mean velocity from the maximum velocity (MXV). Figueras et al. [15] estimated volume flow by single gate Doppler measurements as well. In contrast to Boote [14], Figueras et al. [15] calculated the mean velocity of blood volume flow from the intensity-weighted mean velocity (IWMV) method. Some good experimental results were obtained from their research. Based on the single gate, the narrower the sampling gate, the lower the signal-to-noise ratio (SNR), but the higher the velocity resolution. Correspondingly, the wider the sampling gate, the higher the SNR, the lower the resolution. However, these algorithms based on a single gate do not provide information on the local velocity distribution in the flow cross-section.
This paper proposes a robust ultrasound blood volume flow estimation based on multigate (RMG). The vessel is divided into a plurality of sub-sampling gates. The blood flow velocity in each sub-sampling gate can be accurately estimated as well as the blood volume flow. Therefore, the blood volume flow in each sub-sampling gate can be accumulated to obtain the entire blood volume flow. The multigate method can obtain different blood flow velocities in different location, which provides detail information on the local velocity distribution. When the blood flows in the blood vessel, the blood near the tube wall is affected by viscous friction, the SNR is low so the mean velocity estimation is affected and can may be not obtained accurately. In order to overcome the inaccuracy of velocity, a double iterative algorithm for estimating the mean velocity (DIV) is proposed. The first iteration eliminates the meaningless noise and obtains the accurate signal region. In the obtained signal region, the second iteration integral is made to obtain more accurate estimation of the blood flow velocity. In low SNR conditions, the mean velocity can be estimated accurately using DIV algorithm. In the paper, the RMG method is validated in Doppler phantom and imitation blood flow control system and experimental results are also given. The steadiness in mean velocity estimation for different SNR levels suggests that the DIV algorithm is robust and with low sensitivity to SNR. Through in vitro and in vivo experiments, it can be proved that RMG algorithm has better robustness compared with the existing algorithms, which has important value and significance in clinical diagnostic.
This paper is organized as follows. In Section 2, a detailed description of the algorithm for double iterative flow velocity estimation and multigate blood volume flow estimation is given. Section 3 presents the experimental results and comparisons with the two existing algorithms. Finally, conclusions are drawn in Section 4.
2. Algorithm Description
2.1 Multigate Blood Volume Flow Estimation
In this paper, a robust ultrasound blood volume flow estimation method based on multigate (RMG) is proposed. Multigate technology divides the whole sampling gate into multiple sub-sampling gates. The signal data in each sub-sampling gate is processed respectively. The multigate method can obtain accurate blood flow velocities in each location and provide detail information on the local velocity distribution. Then on the basis of the integral of the vessel cross-sectional area and velocity, the whole volume flow estimation can be obtained accurately.
Assuming that the blood vessel is a circular cross section, as shown in Fig. 1, the entire vessel is divided into even number of sub-sampling gates. The number of gates is N+N' equal to the number of Doppler gates and each sub gate has a size of r1. The signal data in each sampling gate is processed separately and the volume flow at each sampling gate can be estimated.
The blood volume flow Vol(t) at a certain time t through the vessel cross-section is defined as:
(a) The vessel diameter is divided to N+N' number of sub-sampling gates equal to the number of Doppler gates and each gate has a size of r 1. (b) Schematic diagram, assuming that the vessel geometry looks like a circular pipe with a radius of R.
where R is the radius of blood vessel cross-section and V(r, t) is the mean velocity at time t when radius is r.
The discrete form of Eq. (1) at time tj will be
where r1 is the sub gate size, k is the gate index, [TeX:] $$k=1,2, \ldots, N \text { and } V_{k}\left(t_{j}\right)$$ is the velocity estimated in each sampling gate. In Eq. (2) we assume the velocity pattern is symmetric to the vessel diameter and we only count the gate from 1 to N.
In reality, the velocity pattern is asymmetric to the vessel diameter and the vessel diameter is divided to N+N' number of sub gates equal to the number of Doppler gates, as shown in Fig. 1. We can calculate the blood volume flow along the whole vessel diameter, therefore, at time tj , the above equation can be written as
where [TeX:] $$k=1,2, \ldots, N \text { and } k^{\prime}=1,2, \ldots, N^{\prime}.$$
The volume flow estimation by Eq. (3) is just in consideration that the whole sampling gate is divided into even number of sub-sampling gates. Similarly when there are odd number of sub sampling gates N+N'+1, the flow can be estimated as:
Because the ultrasonic scanning beam is not perpendicular to the flow direction, as shown in Fig. 2, so the actual volume flow is a component of estimated volume flow Vol(t) . The volume flow from Eqs. (3) and (4) should be corrected by a factor of [TeX:] $$\sin (\theta),$$ so the actual volume flow [TeX:] $$\overline{V o l}(t)$$ is:
where defines the angle between the ultrasound beam and the flow direction.
The geometry of the ultrasound scanning beam and flow direction.
The heart rate is known to be [TeX:] $$f_{H R}, T_{i-1}$$ represents the end moment of the last cardiac diastolic, and Tj represents the end moment of this cardiac diastolic, then the blood volume flow per minute [TeX:] $$\widehat{V o l}$$ can be estimated as:
The details of multigate blood volume flow estimation algorithm are given by Algorithm 1.
Multigate blood volume flow estimation algorithm
2.2 Double Iterative Flow Velocity Estimation (DIV)
The algorithm proposed for mean velocity estimation is based on the spectral upper and below profile. First find the upper and below profile of the spectral profile image, and then integrate from the upper profile to the below profile. In this paper, a double iterative profile detection algorithm based on energy accumulation function is proposed.
Defining S(n) as an input signal, thus the energy accumulation function P(n) is defined as follows:
where N is the signal length.
The slope of the connecting line between the front and rear of P(n) can be calculated as:
where N is the signal length, P(1) and P(N) are the front and rear of P(n), respectively, S(i) is the input signal, and [TeX:] $$\overline{S}$$ is the mean of signal S(i).
Through simple calculation, it can be known that the slope of the connecting line is similar to the mean [TeX:] $$\overline{S}$$ of signal S(i). Because of the existence of the signal peaks, there must be an evolution process of P(n) from less than [TeX:] $$\overline{S}$$ to more than [TeX:] $$\overline{S}.$$ And in this process, there are three important parameters to be noticed: the negative maximum distance from the connecting line, the positive maximum distance from the connecting line, and the intersecting position.
To get the above three important parameters, Lagrange mean value theorem [16,17] is introduced, that if a function f(x) is continuous on the closed interval [A, B], and is differentiable on the open interval (A, B), therefore there is at least a point on (A, B) to satisfy the following theorem:
In fact, the point of tangency is also the farthest point from the line AB on the curve f(x), namely the point of the positive maximum distance or the negative maximum distance. The energy accumulation curve P(n) satisfies the hypothesis of Lagrange mean value theorem and also in certain SNR conditions, there must be the point of negative maximum distance [TeX:] $$\xi_{m i n},$$ the point of positive maximum distance [TeX:] $$\xi_{m a x},$$ and the intersecting point [TeX:] $$\xi_{\text {cross }}.$$
Rotating the coordinate system of P(n), the rotation angle is defined as [TeX:] $$\theta=\arctan \left(\frac{P(N)-P(1)}{N-1}\right)$$, thus under the new coordinate system, the energy accumulation curve [TeX:] $$\widehat{P}(n)$$ is:
Thus we obtain
Since P(n) is a monotonic increasing function and [TeX:] $$\theta \in(0, \pi / 2),$$ by controlling the order of traversing P(n), it can meet the requirement of [TeX:] $$n \geq \xi_{\min } \text { and } P(n) \geq P\left(\xi_{\min }\right).$$ Therefore, the above complex trigonometric functions can be simplified as follows:
Similarly, the solution formula of [TeX:] $$\xi_{\min } \text { and } \xi_{\text {cross}}$$ can be derived:
[TeX:] $$\frac{P(N)-P(1)}{N-1}$$ is a constant and there is no need to repeat calculating it. The above-mentioned method can avoid the calculation of trigonometric functions brought by coordinate rotation, and it is very meaningful to reduce the computational complexity. Not only that, there is the characteristics of [TeX:] $$\xi_{\min } \leq \xi_{\text {cross }} \leq\xi_{\max },$$ so in real operations, the above three unknown parameters can be got by traversing P(n) only once. In other words, [TeX:] $$\xi_{\max } \text { and } \xi_{\min }$$ are the upper profile and the below profile.
The above process is simply referred to as one iterative algorithm. Because the computational complexity is only O(2N), it is very suitable for real-time systems.
In order to further improve the performance of the algorithm to obtain the details of the profile, this paper designs a double iterative profile detection algorithm. After obtaining the result of one iterative algorithm, we will truncate the P(N) , and then the truncated signal [TeX:] $$\tilde{P}(N)$$ can be obtained:
where Shipt is a constant value of experience. In the strong physical meaning condition (such as blood flow velocity limit), Shipt can be a relatively small value, then one more iteration can obtain more accurate value. On the contrary, in the less physical meaning condition, it requires multiple iterations to find the optimal solution of Shipt. By truncating the signal P(n), the meaningless signal is eliminated, which can improve the SNR of the signal P(n). The truncated signal
[TeX:] $$\tilde{P}(N)$$ performs the same processing as one iterative algorithm, and further obtains the local optimal solution, thus improves the detailed resolution of the profile.
Generally, intensity-weighted integral is used to calculate mean velocity. Thus the mean velocity [TeX:] $$V_{k}(t) \text { of the } k^{\text {th }}$$ gate at time t could be given by
where [TeX:] $$P_{t}^{k}(x)$$ is Doppler power spectrum of the [TeX:] $$k^{t h}$$ gate at time t, and V(x) is the velocity of the [TeX:] $$k^{t h}$$ gate.
However, intensity-weighted integral is a global algorithm and can ignore local details. To overcome the inaccuracy, this paper proposes an accurate flow velocity estimation algorithm based on the double iterative profile detection. The double iterative algorithm changes the range from the global to the local, and eliminates the unnecessary information that interferes with the signal, so the accuracy will increase. In the estimation of mean velocity, the calculation range is calculated from the upper profile position to the below profile position.
Thus the mean velocity [TeX:] $$V_{k}(t) \text { of the } k^{\text {th }}$$ gate at time t could be given by
where [TeX:] $$P_{t}^{k}(x)$$ is Doppler power spectrum of the [TeX:] $$k^{t h}$$ gate at time t, [TeX:] $$S_{u p p e r}(t) \text { and } S_{\text {below}}(t)$$ are the upper profile position and the below profile position at time t, respectively, and V(x) is the velocity of the [TeX:] $$k^{t h}$$ gate.
The velocity pattern V(x) at time t.
As shown in Fig. 3, the power spectrum within the upper and the below profile in each sampling gate will be used to calculate the mean velocity of Eq. (17).
The details of double iterative flow velocity estimation algorithm are given by Algorithm 2.
Double iterative flow velocity estimation algorithm
3. Experimental Results
The proposed DIV algorithm is compared with the two existing algorithms, IWMV and MXV. The three algorithms for mean velocity estimation are tested in vivo recordings. And also, flow simulations are performed and the algorithms are tested on them. In this section, the accuracy of mean velocity estimation is evaluated and the robustness of the proposed algorithm is presented.
The volume flow estimation based on the IWMV and MXV of single gate is presented, referred as intensity-weighted mean flow (IWMF) and maximum flow (MXF). MXF is defined as maximum flow. In the proposed RMG algorithm, the blood vessel is split into 10 sub-sampling gates. According to the double iterative velocity estimation formula, the mean velocity of the 10 sampling gates can be calculated respectively, as well as the volume flow. Then the volume flow of 10 sampling gates can be accumulated to get the whole vessel volume flow.
The in vitro tests have been carried out through the experimental set-up, KS205D-1 type Doppler phantom and imitation blood flow control system developed by the Institute of Acoustics, Chinese Academy of Sciences. This set-up can measure the blood volume flow. With this set-up, the performance of the algorithms can be verified.
The measurements performed with the same pump setting were compared for different methods. The error is the ratio between the absolute error caused by the measurement and the actual value of the volume flow. As shown in Table 1, under different flow settings, estimated volume flow values and mean errors by the RMG method and the two comparison methods (IWMF and MXF) are presented.
As shown in Table 1, at different flow rates, the error of the proposed RMG algorithm is significantly lower than the MXF and IWMF. The RMG algorithm can more accurately measure the volume flow under different flow settings. The mean error of the proposed RMG algorithm is 4.81%, which is much smaller than the mean error obtained by the traditional single sampling gate algorithms MXF and IWMF. Due to the steadiness and robustness of the RMG algorithm, it can be used in clinical practice to estimate volume flow and detect cardiovascular diseases.
Experiments to test the performance of the DIV algorithm were done for several different SNR levels. The data we deal with in this sectuib comes from the radio frequency (RF) carotid blood flow signals collected by ultrasound instruments. On the basis of the original spectral profile image, we add the Gaussian white noise with a constant mean of 0 and variance of 0.01, 0.02, 0.03, 0.04, and 0.05, respectively, to test the performance of the three algorithms in different SNR levels. The ideal curve of the experiment is the mean velocity curve drawn by an experienced clinician according to the spectral profile image.
Estimated volume flow values and mean errors by the proposed method and the two comparison methods
Fig. 4 shows the experimental results of mean velocity estimation in the spectral profile images of six different SNR conditions. It is found that in the lower SNR spectral profile images, as shown in Fig. 4, the IWMV and MXV have larger errors. With the increase of noise, the mean velocity curves obtained by IWMV and MXV show serious jitter and it leads to inaccurate calculation of the mean velocity at the local details. Therefore, in the results from IWMV and MXV in Fig. 4, the mean velocity estimation is failing in the case of greater noise. In the same condition of low SNR, the proposed DIV algorithm works well.
where M is the actual number of estimated mean velocity points, N is the number of standard mean velocity points, is the scaling constant (usually 1/9), and d(i) is the distance between estimated mean velocity point and standard mean velocity.
Experiments on 10 sets of data are carried out to obtain the mean values and standard deviations of PFOM values under different noise conditions. We add the Gaussian white noise with a constant mean of 0 and variance of 0.01, 0.02, 0.03, 0.04, and 0.05, respectively, to each set of original data. Fig. 5 gives the results of PFOM mean for the three algorithms. The PFOM mean of DIV is significantly higher than the IWMV and MXV. The PFOM mean of DIV is closer to 1 and it shows that the mean velocity curve obtained by DIV is closer to the ideal curve drawn by the experienced clinician. Moreover, the line chart shows that the PFOM mean of the algorithm DIV is more stable and its jitter is smaller in the case of different noises.
Comparison of the mean velocity estimation using IWMV, MXV, and DIV algorithms and the ideal curve drawn by an experienced clinician, on carotid artery spectral profile images, in different noise conditions: (a) the original spectral profile image, (b)–(f) are the spectral profile images added Gaussian white noise with mean of 0 and variance of 0.01, 0.02, 0.03, 0.04, and 0.05, respectively.
Fairly low standard deviations of DIV algorithm in different noise conditions are found in Table 2. Less than 6% mean standard deviation by using DIV algorithm to estimate mean velocity is observed. Therefore, it is steady when Gaussian white noise is added with mean of 0 and variance range of 0.1–0.5, which indicates that the mean velocity estimated by DIV has very low sensitivity to SNR. Mean velocity estimated using DIV in vivo data can be observed to show better performance compared with IWMV and MXV at the same SNR level. These results point toward robust mean velocity estimation of DIV algorithm with low sensitivity to SNR. When the noise level is very high, the algorithm MXV fails completely and the PFOM mean of each set of data is very small, as shown in Fig. 5, so the standard deviation is very small. Therefore, there will be such a case that the standard deviation is small in Table 2.
PFOM mean of DIV, IWMV, and MXV algorithms by varying variance of Gaussian white noise from 0.01 to 0.05 when mean of Gaussian white noise is 0.
4. Conclusion
In this paper a robust ultrasound blood volume flow estimation algorithm based on multigate (RMG) and an accurate double iterative algorithm for estimating the mean velocity (DIV) have been presented. Based on the multigate method, the velocity and volume flow of each sub sampling gate can be estimated, which provides clinicians with detail information on the local velocity distribution. Experiments to test the performance of the DIV algorithm were done for different SNR levels. Low errors in mean velocity estimation are observed. Mean standard deviation of the mean velocity estimation is 5.34% and is significantly low. The stability in PFOM calculation for different SNR levels suggest that the DIV algorithm is robust and with low sensitivity to SNR. The RMG algorithm is tested in a custom-designed experimental setup, Doppler phantom and imitation blood flow control system. Less than 5% mean error suggests that the RMG algorithm provides better performance compared with the existing volume flow estimation algorithms. Estimation of accurate blood volume flow in ultrasound Doppler blood flow spectrograms is of important guiding significance for monitoring cardiovascular diseases. Therefore, the stability and robustness of the proposed algorithm shows that it can be widely applied to real-time blood volume flow estimation in clinical practice, which can be useful in monitoring cardiovascular diseases. In summary, we are optimistic about the potential of this new approach to blood volume flow estimation for use in future medical instruments.
Acknowledgement
The authors would like to thank the engineers at the Saset (Chengdu) Inc. who contributed to the experimental works done in Saset’s commercial ultrasound machine (iNSIGHT 23R).