A robust local pulse wave imaging method based on digital image processing techniques

1.   Introduction

Currently, cardiovascular disease (CVD) is the major cause of mortality ^[1]. According to the China Cardiovascular Disease Report 2018, the prevalence and mortality of cardiovascular disease in China are still rising, with a projected 290 million people with current cardiovascular disease ^[2]. Early detection and risk assessment of CVD risk factors can help to reduce the overall CVD-related morbidity and mortality ^[3,4]. Atherosclerosis is widely associated with cardiovascular risk and provides a powerful signal for assessing cardiovascular risk ^[5].

Pulse wave velocity (PWV) has been demonstrated to correlate with age, coronary artery disease, and hypertension, making it one of the most dependable markers for assessing vascular sclerosis ^[6,7,8]. However, the majority of PWV utilized in this sort of research are brachial-ankle, cervical-femoral, and femoral-ankle PWV ^[9,10]. Their calculation of the distance between two points relies on the assumption that the vessels in the body are aligned in a straight line, which introduces measurement inaccuracy. Moreover, this measurement method reflects the average PWV between two detection points, with atherosclerosis and abdominal aortic aneurysm limitations ^[11].

Local PWV measurement of arteries can effectively compensate for the deficiencies of the above PWV measurement methods ^[12]. Qian and Li et al. ^[13,14] published studies on the correlation between local PWV of carotid arteries and sex, age, and intima-media thickness. Their experiments were performed with VINNO M80 ultrasound equipment, which integrates a technique based on a fast comb-tooth-shaped focus beam transmitting and receiving sequence local PWV detection. However, this technique also has limitations regarding arterial tortuosity ^[15]. Recently, Yang et al. ^[16] examined the association between diabetes mellitus and carotid PWV using the Supersonic Imagine Aixplorer Color doppler ultrasonic diagnostic apparatus from the Supersonic Imagine company. This study further expands the scope of clinical applications of PWV. Speckles in ultrasound imaging affect image quality and can make post-processing difficult. Liu et al. ^[17] proposed a robust detail-preserving anisotropic diffusion filter for speckle reduction in ultrasound images. Tang et al. ^[18] develop a new bilateral filter for speckle reduction in ultrasound images for follicle segmentation and measurement.

In studying local PWV measurement methods, Luo's group at Tsinghua University conducted an early study of the pulse wave imaging (PWI) technique. Moreover, this technique realized the visualization of PWV propagation. In 2012, they measured the carotid PWV values of eight healthy young people using the PWI method to verify the feasibility of the method ^[19]. In 2014, ^[20] utilized the PWI method for tests on simulated data to study the effect of ultrasound imaging's frame rate, number of scan lines, and picture width on the precision of PWV value assessment, which provides an essential reference for the optimization of local PWV measurement parameters.

Local PWV measurement techniques can be classified into original radio frequency (RF) signal-based measurements and ultrasound image-based data measurements ^[21]. The first method is similar to the PWI method mentioned above, but this method leads to the measurement of local PWV only by specific ultrasound machines. Negoita et al. ^[22,23] applied ultrasound images to measure utilizing the ln(D)U-loop method, which belongs to the second method, and its measurement accuracy has some differences compared to the first method but has the advantage that it can be independent of the ultrasound imaging system and can be easily promoted and widely used. Rabben et al. ^[24] were the first to implement the QA-loop method to the human and canine carotid by utilizing ultrasonography to determine blood flow volume and cross-sectional area of blood arteries and comparing it to the foot-to-foot method. In 2010, Feng et al. ^[25] used the ln(D)U-loop method to measure the diameter of arteries and the variation in blood flow velocity by measuring local PWV of the carotid artery in humans and the ascending aorta in dogs. Negoita et al. ^[22,23] first measured the PWV of the ascending aorta in 13 healthy humans using this method based on ultrasound images, but their proposed method still has problems such as high image quality requirements when extracting the diameter waveform and a narrow range of experimental data. In this paper, we attempt to address the problem that the ln(D)U-loop method for measuring the PWV of the ascending aorta requires an excessive amount of ultrasound image quality and manual intervention. We improve the diameter waveform extraction method by combining the gradient calculation and the center of gravity method ^[26]. At the same time, the linear growth part of cardiac systole was found using piecewise linear regression in this paper so that the average value of the slope obtained from its calculation can more accurately calculate PWV values.

2.   Theoretical analysis

2.1.   ln(D)U-loop method

The ln(D)U-loop method was proposed by Feng and Khir ^[25]. Wave speed is dependent on the distensibility of the tube wall, ${D}_{s}$, and the density of the fluid, $\rho$

where ${D}_{s}$ = (1/A) (dA/dP), dA is the change in a circular cross-sectional area and dP is the change in pressure. Therefore, substituting the formula for the area of a circle in Eq (1) gives

where D is the diameter. The well-established pressure equation introduced by Kries et al. ^[27] was written as

Substituting Eq (3) into the original PWV equation gives

If we consider substituting dD/D = dlnD into Eq (4), the PWV representation in the ln(D)U-loop method gives

2.2.   Proposed method

The flow of the ln(D)U-loop method is shown in Figure 1.

Figure 1.  The flow of the ln(D)U-loop method.

DownLoad: Full-Size Img PowerPoint

Step 1: The peak of R-waves in the electrocardiogram (ECG) is used to separate individual heartbeats.

Step 2: Extract the diameter waveform from the M-mode ultrasound image of the ascending aorta over the duration of one heartbeat.

Step 3: Extract the velocity waveform from the PW-mode ultrasound image of the ascending aorta over the duration of one heartbeat.

Step 4: Determine the systolic foot (SF) of both the extracted diameter waveform and velocity waveform as the starting point of the cardiac systole and calculate the linear part of the loop.

2.2.1.   Diameter waveform extraction

Using the ECG in the M-mode ultrasound image, the individual heartbeat is divided according to the R-wave peak. Figure 2 illustrates the workflow for obtaining the diameter waveform. The gradient filtering of abnormal values is as follows:

Figure 2.  Flow chart of extracting the diameter waveform from the M-mode ultrasound image.

DownLoad: Full-Size Img PowerPoint

where P is the set of outlier points after filtering, and G is the threshold of the gradient. $\left|\text{g(i)}\right|$ is the absolute value of the i-th gradient.

Figure 3(a) depicts the result of superimposing adjacent cardiac cycles in the M-mode ultrasound image and filtering away the majority of the noise using the thresholding technique. However, the top and bottom vessel walls that are roughly found still contain abnormal points. Figure 3(b) depicts the result of using gradient filtering to filter out the abnormal points. Figure 3(c) depicts the second search results in the top and bottom vessel walls using the center of gravity method for the thresholded M-mode ultrasound images. Figure 3(d) shows the image after processing by Eq (8). $wal{l}_{T}$ in Eq (7) is the top and bottom vessel walls found by the thresholding technique, $wal{l}_{G}$ is the top and bottom vessel walls found by the center of gravity method, and $\alpha$ is the adjustment scale factor, which takes a value between 0 and 1.

Figure 3.  The ascending aorta diameter waveform extraction process. (a) Anterior and posterior wall of blood vessel obtained by threshold method (b) The result of the gradient of (a) after filtering out outliers (c) The center of gravity method to find the position of the blood vessel wall (d) Results of (b) and (c) using coefficient stacking (e) After spline smoothing, the diameter of the blood vessel subtracted from the anterior and posterior walls, the red dot indicates the constricted foot of the found diameter D waveform.

DownLoad: Full-Size Img PowerPoint

Finally, interpolation smoothing is performed using a smoothing spline for the top and bottom vessel walls in Figure 3(d). Subtracting the positions of the two walls generates the diameter waveform depicted in Figure 3(e). The x-axis indicates the cardiac cycle time calculated from the heart rate, and the y-axis indicates the diameter converted from scales and pixels.

The linear regression method for SF in this paper is similar to that of Negoita ^[22]. As shown in Figure 3(e), the red point represents the starting point of the obtained diameter waveform.

In addition, a set of experiments achieves the finding of anterior and posterior walls by the thresholding method mentioned above. In Figure 4(a)–(f) are images with a good visual effect of thresholding between 75 and 85 with an adjacent interval of 2. The arrows in the images mark the points where there are still some abnormalities after the thresholding segmentation. Figure 4(g) shows the positions of the anterior and posterior walls of the blood vessels detected by our method. It can be observed that without the overlay of adjacent images, achieving a good result using only the constant threshold adjustment is not easy.

Figure 4.  Comparison of different thresholds with the methods used in this paper to detect the anterior and posterior walls of blood vessels. (a)–(f) are images with a good visual effect of thresholding between 75 and 85 with an adjacent interval of 2. (g) shows the positions of the anterior and posterior walls of the blood vessels detected by our method.

DownLoad: Full-Size Img PowerPoint

2.2.2.   Velocity waveform extraction

According to pulsed wave Doppler image characteristics of the ascending aorta, its pixels present higher brightness at both ends of the cardiac cycle. A threshold of 4 times the middle part was used for the left 1/8 pixels of the pulsed Doppler image for one cardiac cycle, and a threshold of 3 times the middle part was used for the right 2/5 pixels of the image. This parameter is adjusted according to the characteristics of the image and has better results in this experiment. The user can modify this parameter appropriately according to the visual effect. After the improved thresholding to obtain the maximum envelope, Figure 5(a) shows the image obtained by filtering the anomalous outliers using the gradient method. Finally, as shown in Figure 5(b), a blood flow velocity curve of the cardiac cycle over time is obtained by spline smoothing. The coordinate conversion method is obtained by transforming the scale concerning pixels. Similarly, as shown in the red dots in Figure 5(b), a linear regression method was used to find the SF location of the blood flow velocity waveform.

Figure 5.  The ascending aorta velocity waveform extraction process. (a) The result of gradient filtering on the maximum envelope obtained by the improved threshold method (b) Velocity waveform obtained by applying smoothing splines to the envelope of (a), the red dots indicate the contraction legs of the velocity U waveform.

DownLoad: Full-Size Img PowerPoint

2.2.3.   Segmented linear regression for PWV detection

In the ln(D)U-loop method theory, PWV needs to calculate the linear growth part in the absence of reflections. In this paper, we propose to use a segmented linear regression method for fitting. Firstly, we select 20 data points as the first window, find a line with maximum linear regression r2 from the first 1 to 10 points, denoted by ${\mathrm{y}}_{1}$ = ${\mathrm{a}}_{1}$x+${\mathrm{b}}_{1}$, and then find a line with maximum linear regression r2 from 10 to 20 points, denoted by ${\mathrm{y}}_{2}$ = ${\mathrm{a}}_{2}$x+${\mathrm{b}}_{2}$.

Finally, as shown in Eq (8), the average of the two straight lines is the linear growth part of the systolic period to be found, using the average of the slopes to calculate PWV.

The PWV of all combinations was sorted and filtered out of the front and back 25% of possible outliers, and the mean (± standard deviation) was used to calculate the PWV of volunteers. We employed random sample consensus (RANSAC) for the absence of reflections fitting the loops and compared it to the original single-linear linear regression to analyze the performance of the piecewise linear regression approach suggested in this paper.

The original single linear regression method is to perform a linear regression after the SF point until the r² of the linear regression is less than some set value. The disadvantage of this method is that a threshold value of r² needs to be given, which affects the accuracy of the linear fit and the calculation of PWV. Especially in the case of poorly extracted velocity and diameter waveforms, the anomalies can easily lead to large deviations in the fit. RANSAC is fitted by given mathematical expressions. The division of the cardiac cycle is not just systolic and dilated, and a more detailed segmentation of the cycle is beneficial to the experiment.

3.   Experimental setting

3.1.   Data collection

In this study, a total of 12 healthy young volunteers (age 32.5 ± 5.6, mean ± standard deviation) and 17 elderly patients (age 60.5 ± 9.2) with cerebral infarction were recruited and scanned with an ultrasound system, GE Vivid E95 (GE Healthcare: Vingmed Ultrasound, Horten, Norway) ultrasound device. M-mode ultrasound and PW-mode ultrasound images were acquired from the ascending aorta in all volunteers at the same location with the same scan speed. Both ultrasound modes had 3-lead ECG recordings. Each image had four to five R-wave peaks, and two to three images of different modes were acquired for each volunteer.

3.2.   Evaluating performance

Each volunteer had 16 sets of measured PWV. After sorting and eliminating the 16 sets of measured PWV, the mean was calculated and reported as the mean ± standard deviation ($\overline{PWV}$ ± SD). Usually, the true PWV value is not available in somatic measurements, so it is not possible to directly evaluate the accuracy of the measurement by the improved method in this paper. The statistical difference between the two sets of measurements was measured using the Wilcoxon rank-sum test, and a P value less than 0.01 was defined in this study as a difference that was considered statistically significant.

To further evaluate the method proposed in this paper, the QA-loop method of measuring PWV was used for comparison. It matches blood flow Q and cross-sectional area A using the same rules as the ln(D)U-loop method, and the postprocessing is also consistent with the ln(D)U-loop method. Blood flow Q was calculated using the method in Holdsworth ^[29].

4.   Results

The effectiveness of the experimental results is evaluated by the ability to accurately fit the data points in the rising period as a judgment. The fitting results of RANSAC and the segmented linear regression proposed in this paper in which different results appear are shown in the blue line in Figure 6. Figure 6(a) shows the fitting results of RANSAC, and Figure 6(b) shows the fitting results of segmented linear regression. Due to the offset of SF point finding of some test objects, this part of the anomaly can lead to bias in the fitting effect of RANSAC algorithm. From the comparison of Figure 6(a) and (b), the segmented linear regression method proposed in this paper is relatively non-sensitive to some anomalies and can achieve a good fitting effect. Note that the first 10 data points in Figure 6(b) are fitted with the y₁, the last 10 data points are fitted with the y₂, and the final fitted line for y₁ and y₂ is the blue line y.

Figure 6.  RANSAC and piecewise linear regression are two methods of fitting effect comparison. (a) Using RANSAC to fit the results of the linear part (b) Piecewise linear regression method to fit the result of the linear part.

DownLoad: Full-Size Img PowerPoint

The blue line in Figure 7 shows the linear part of the automatically calculated loop. It compares the original single linear regression method with the two linear regression averaging methods proposed in this paper. Figure 7(a) shows that the single linear regression method cannot fit the data points in the rising period better. The improved linear regression method shown in Figure 7(b) can provide better coverage of the data points in the rising period.

Figure 7.  Comparison of the results of the original linear regression and the improved linear regression method. (a) Linear part plot of the original method with R₂ of 0.97 (b) The linear part graph after the improved linear regression in this paper.

DownLoad: Full-Size Img PowerPoint

The PWV measurement results of each volunteer by the QA-loop method and ln(D)U-loop method are shown in Table 1. The first 12 groups represent the findings of a healthy group with a mean PWV = 6.81 ± 1.73 $m{s}^{-1}$, and the last 17 groups represent the results of elderly patients with PWV = 12.56 ± 3.47 $m{s}^{-1}$. There is a considerable disparity between the PWV measurements of the two groups. The ln(D)U-loop and QA-loop methods have P-values of 0.00029 and 0.001, respectively, indicating that the PWV calculated by both methods can effectively screen cerebral infarction patients and normal populations.

Table 1.  PWV measurement results of each volunteer by the QA-loop method and ln(D)U-loop method.

Volunteer No.	Cerebral Y/N	$PW{V}_{QA}$/$m{s}^{-1}$	$PW{V}_{ln\left(D\right)U}$/$m{s}^{-1}$
1	NO	9.08 ± 0.59	8.09 ± 0.87
2	NO	5.25 ± 0.52	4.52 ± 0.38
3	NO	7.45 ± 1.38	7.62 ± 1.54
4	NO	6.36 ± 0.59	6.38 ± 0.80
5	NO	9.01 ± 0.24	9.07 ± 0.23
6	NO	6.32 ± 0.33	6.67 ± 0.61
7	NO	6.87 ± 0.79	7.51 ± 0.67
8	NO	8.01 ± 1.23	6.95 ± 0.42
9	NO	8.80 ± 0.38	8.04 ± 0.37
10 11 12	NO NO NO	2.90 ± 0.50 5.67 ± 0.25 7.14 ± 0.14	3.31 ± 0.81 6.41 ± 0.10 7.15 ± 0.04
13	YES	9.01 ± 0.55	9.02 ± 1.01
14	YES	18.37 ± 1.09	18.38 ± 0.93
15	YES	13.62 ± 1.06	14.35 ± 1.94
16	YES	17.79 ± 3.32	17.28 ± 2.76
17	YES	9.54 ± 1.59	8.79 ± 0.84
18	YES	8.89 ± 1.86	9.53 ± 1.07
19	YES	10.58 ± 1.05	10.38 ± 0.96
20	YES	8.40 ± 1.10	9.54 ± 1.74
21	YES	11.23 ± 2.71	11.46 ± 3.79
22	YES	8.60 ± 1.40	8.77 ± 1.48
23 24 25 26 27 28 29	YES YES YES YES YES YES YES	13.81 ± 2.15 9.91 ± 0.09 16.03 ± 0.07 15.75 ± 0.36 18.49 ± 0.02 9.94 ± 0.03 11.32 ± 1.06	13.76 ± 3.16 10.16 ± 0.16 16.50 ± 0.07 15.85 ± 0.73 18.01 ± 0.03 10.79 ± 0.04 11.01 ± 0.88

---

 | Show Table

DownLoad: CSV

An analysis of the agreement (Bland-Altman analysis) between the QA-loop method and ln(d)U-loop method. is shown in Figure 8. The mean deviation of the measured PWV values is 0.07 m/s, and the standard deviation of the deviation is 1.18 m/s.

Figure 8.  Agreement analysis of the measured PWV values between the QA-loop method and the ln(d)U-loop method.

DownLoad: Full-Size Img PowerPoint

5.   Discussions and conclusions

Numerous clinical investigations have established that noninvasive local PWV measurement methods accurately represent arterial stiffness information ^{[13,14,15,16]}. This study focuses on improving the utilization of ultrasound image data required for the ln(D)U-loop method and reducing manual intervention. The results show that for the rising point data fitting (Figures 6 and 7), neither the single linear fit nor the RANSAC method can fit the data better. This shows that the strategy provided in this paper is efficient.

In addition, although there is some deviation between the PWV measured by the two methods, the standard deviation of the deviation is relatively small. Considering the mean difference in PWV between the healthy and infarct groups of 5.75 m/s and from the Wilcoxon rank-sum test, both methods can distinguish well between the two data groups.

It is important to note that the differences reflected by PWV measurement include multiple factors, such as age and blood pressure. Many other physiological factors affecting PWV are not considered. These are very relevant clinical issues but have higher requirements for sample size, which are left for follow-up studies.

In this study, an improved ln(D)U-loop method is proposed to measure PWV, and the poor-quality image data of M-mode ultrasound can be processed using a diameter waveform extraction method that combines thresholding, gradient filtering, and the center of gravity method. The piecewise linear regression method also calculates PWV to reduce errors. This study increases the possibility that the ln(D)U-loop method will be extensively employed in the clinical evaluation of local PWV, which has the potential to become a versatile diagnostic adjunct for cardiovascular diseases such as hypertension.

Acknowledgments

We thank the editors and anonymous reviewers for their helpful comments for improving our work. The study is funded in part by a research grant from the Science and Technology Council of Sichuan Province (2021YJ0248), and Natural Science Foundation of Sichuan Province (2022NSFSC0833).

Conflict of interest

All authors declare no conflicts of interest in this paper.