Keywords

1 Introduction

Pulmonary embolism (PE) is a blockage in the lungs that pulmonary artery (PA) is obstructed by emboli, leading to pulmonary circulatory disturbance with low diagnosis rate and high mortality. The emboli include thrombosis, fat, amniotic fluid, air and tumor, among which thrombosis is the most common. PE usually results from a blood clot in the leg where the thromboses fall off and move with blood circulation into PA. Therefore, the pulmonary embolism only occurs in the pulmonary artery and its branches. Symptoms of a PE include dyspnea, chest pain and hemoptysis. PE can be divided into acute pulmonary embolism (APE) and chronic pulmonary embolism (CPE) according to the size and location of emboli. APE is the most common type of PE. According to the American Heart Association (AHA) 2015 data, PE is the third leading cause of cardiovascular death. It is difficult to diagnose PE definitively from other causes of symptoms for lack of clinical specific manifestation. The risk assessment is classified according to the pulmonary artery severity index, including the percentage of emboli volume to PA volume, the location of the emboli in the pulmonary artery branches, the spatial relationship between the emboli and the PA, and the type of emboli. They are all deemed as standard risk assessment indicators. We need computers to assist people to predict the risk of individual illness or death in diagnosing PE [1]. Since the PE only occurs in the pulmonary artery and its branches, accurate segmentation of the PA is the primary task.

Spiral CT pulmonary angiography (CTPA) is the golden standard for the PE diagnosis. CTPA images can clearly show the location, shape, size of the thromboses in the PA subsegment branches or above [2]. The sensitivity and specificity of the diagnosis are accurate, and it has many merits including non-invasive nature, fast scanning speed and powerful image processing ability. These advantages are instructive for clinical diagnosis and treatment plan, especially suitable for APE patients in emergency.

Medical image segmentation refers to extracting regions of interest according to the certain features or similarity of feature sets in medical images, and the images are divided into several non-overlapping regions with certain consistency. Traditional methods comprise threshold-based methods (Otsu method, etc.), region-based methods (Region Growth method, RG, etc.), boundary-based methods (Canny operator method, etc.). Kanas et al. proposed a random walk algorithm based on local intensity differences to segment the explicit boundaries of brain tumor [3]. Lesage D. et al. proposed a Bayesian random tracking algorithm based on particle filter for segmentation of coronary artery trees from cardiac CTA data. Using kernel density estimation to learn the entire Bayesian model, the results turned out to be a high robustness [4]. Zhang et al. proposed a piece-wise linear transformation for enhancing brain vessel boundaries and then segmenting blood vessels. The blood vessel was fused with weighted method to gain the entire vessels [5]. Zhou et al. presented a three-dimensional voxel clustering with multi-stage adaption method to segment pulmonary vessels based on expectation-maximization (EM) analysis. After analyzing connectivity of region on the segmentation results, blood vessels were tracked and reconstructed [6]. Some new methods, such as fuzzy mathematics-based methods (Fuzzy C-means algorithm, FCM) and neural network-based methods (Convolutional Neural Network, CNN), are put forward in later studies. Qin et al. proposed a convolutional-deconvolution depth network model with residual connections for segmenting the lesions of prostate in MR images. Combination of the U-Net and the ResNet network achieved a better result [7]. Xu et al. presented a new stage-wise convolution network for segmenting pulmonary vessels. The network was characterized by the automatic learning of pulmonary vessels in stages [8]. Pulmonary vascular segmentation is often used to detect pulmonary nodules [9], so the researchers did not separate the pulmonary artery and pulmonary veins.

The region growth (RG) method based on region segmentation has been widely applied to blood vessel segmentation, including two-dimensional RG and three-dimensional RG. The method fully considers the connectivity of blood vessels in space. In the study of Region Growth algorithm for pulmonary artery segmentation, scholars have proposed some improved methods. Zhang segmented the pulmonary trunk using a method based on RG and slice marching. The branches were tracked from the pulmonary trunk to branches, however, there existed fracture in branches and roughness in vascular wall after reconstruction [10]. Ebrahimdoost et al. took the pulmonary trunk extracted by the RG as the initial contour of the level set segmentation. They removed adhesions from other tissues and finally the pulmonary artery could be classified to the lung segment level [11]. Flores et al. proposed a 3D region growth plus progressive method on PA segmentation. They introduced boundary-based termination conditions, which could segment the PA to the segment level [12].

This paper proposes an improved three-dimensional region growth approach to segment pulmonary artery in the following three parts. Firstly, to separate vena cava, on the basis of mathematical principle and anatomical features, the vector of the pulmonary trunk in three-dimensional space is obtained. An isolation plane is introduced to remove vena cava. Secondly, in order to reduce error in different doctor’s experience, the seed point is automatically obtained. Thirdly, to solve the precise segmentation to subsegmental branches, an adaptive parameter is introduced. The CT threshold will iterate according to the direction in which the CT value decreases. The approach for PA segmentation is presented in Sect. 2, and the results and conclusion are shown in Sect. 3 and 4, respectively.

2 The Proposed Method

In this paper, we fully take into account the connectivity nature of three-dimensional region growing algorithm. The false positives detected by the existing pulmonary embolism system mainly occur in the pulmonary vein and there is adhesion between the branches of pulmonary artery and pulmonary vein. So, the idea of connected domains is applied to address the problem. And in this part, the proposed 3D RG algorithm is illustrated in detail.

2.1 Remove Vena Cava

The blood from the vena cava flows through the heart into the pulmonary artery and there is little difference in CT value between the coterminous tissues. If the three-dimensional region growth is carried out directly, two regions will be segmented together. That is, the phenomenon of leakage occurs. In order to avoid leakage, we find the junction of pulmonary trunk and vena cava, then remove vena cava. We obtain normal vector of the pulmonary trunk in cross-section plane and add an isolation plane at the junction. Since obtaining the three-dimensional vector is difficult, we start from the two-dimensional vector. In CTPA images, the outline of the pulmonary trunk looks like a Chinese character- ‘人’. The main pulmonary artery (MPA) connected to the vena cava extends upward, the left pulmonary artery (LPA) and the right pulmonary artery (RPA) extend downward. First, find the pulmonary trunk. Then, extract the skeleton of pulmonary trunk to gain the vector. Lastly, calculate the three-dimensional vector at the junction according to the vector projection theorem in mathematics. The detailed algorithm is as follows:

The first step, find the pulmonary trunk in the slice. We select a CTPA slice containing PA region (Fig. 1(a)). The major highlighted region on the image segmented with threshold method is the pulmonary trunk (Fig. 1(b)).

Fig. 1.
figure 1

(a) The original image. (b) The pulmonary trunk in the slice. (c) The skeleton. (d) Geometric diagram of skeleton.

Full size image

The second step, extract the skeleton of pulmonary trunk in the slice. The process of extracting the skeleton is performed in the region of ‘人’. Extraction Skeleton of the region uses a thinning algorithm in morphological principle called Medial Axis Transform (MAT) [13]. In the region R, the set of all elements is A. For any point p in set A, find the closest point on the boundary b. If there exist more than one closest point, it is considered that p is on the medial axis of the region R. So, the set of p is the skeleton of the region. The region R is eroded by the structuring element B, and the skeleton expressions are shown in the Eqs. (1) and (2):

$$ S(A) = \bigcup\limits_{k = 0}^{K} {S_{k} (A)} $$
(1)

And

$$ S_{k} (A) = (A\Theta kB) - (A\Theta kB) \circ B $$
(2)

Where: \( A\Theta kB \) represents eroding \( k \) times to \( A \), and \( k \) is the last time before being eroded into an empty set.

After pulmonary trunk is eroded \( k \) times according to Eqs. (1) and (2), the skeleton is obtained as shown in Fig. 1(c). So, we get the skeleton center point O and the endpoint P in MPA. And then, a vector \( \overrightarrow {OP} \) or \( \overrightarrow {{n_{0} }} \) (Fig. 1(d)) is obtained from point O and P, which is the vector in the slice of the MPA. The process of skeleton extraction is shown in Fig. 1. The original image is shown in Fig. 1(a). The pulmonary trunk in the slice is shown in Fig. 1(b). And the skeleton is shown in Fig. 1(c). The geometric diagram of skeleton is shown in Fig. 1(d).

The third step, get the spatial normal vector of MPA. According to the vector projection theorem, \( \overrightarrow {n} \) is the vector of the three-dimensional space, and \( \overrightarrow {{n_{0} }} \), the projection of \( \overrightarrow {n} \) onto the XOY plane, is the vector of two-dimensional plane. \( \left| {\overrightarrow {{n_{0} }} } \right| \) is the length of \( \overrightarrow {{n_{0} }} \), and \( \left| {\overrightarrow {n} } \right| \) is the length of \( \overrightarrow {n} \). θ is the angle between \( \overrightarrow {{n_{0} }} \) and \( \overrightarrow {n} \) in the space. The theorem expression is shown in the Eq. (3):

$$ \left| {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{{n_{0} }} } \right| = \left| {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{n} } \right|*\cos \theta $$
(3)

The pulmonary trunk is a spatial region at a certain angle from the XOY plane in the three-dimensional space (Fig. 2(a)). Figure 2(b) is the right view of Fig. 2(a). The pulmonary trunk, a short and wide vessel, ascends obliquely, being divided into LPA and RPA [14]. In our study, \( \overrightarrow {n} \) is spatial normal vector of three-dimensional MPA in cross-section plane, XOY plane is horizontal plane, and \( \overrightarrow {{n_{0} }} \) is the vector of the MPA in the XOY plane which has been obtained in second step. θ, between three-dimensional MPA and the XOY plane, is known according to anatomical features. And then vector \( \overrightarrow {OQ} \) (\( \overrightarrow {n} \)), the actual MPA vector, is acquired (Fig. 2(c)). The 3D reconstruction result of pulmonary trunk is shown in Fig. 2(a), the right view of pulmonary trunk is shown in Fig. 2(b), and geometric diagram is shown in Fig. 2(c).

Fig. 2.
figure 2

(a) The 3D pulmonary trunk result. (b) The right view of pulmonary trunk. (c) Geometric diagram of MPA vector.

Full size image

The fourth step, add the isolation plane at the junction. The MPA is approximately 5 cm in length and 3 cm in diameter [14]. The MPA endpoint W can be obtained from the point Q and vector \( \overrightarrow {n} \) (Fig. 3(a)). \( \left| {\overrightarrow {OW} } \right| \) is the length of MPA. It is confirmed that the point W must be in the pulmonary artery region. In space, MPA is similar as cylinder. With the endpoint W being the center and \( \overrightarrow {n} \) being the normal vector, create a small circular plane that is the isolation plane for the three-dimensional region growth. The diameter of the circle is slightly larger than the MPA. This algorithm is suitable for angle between 0°–90°, except 0° and 90°. Figure 3(a) is the geometric diagram of the spatial position of the isolation plane. Figure 3(b) is a slice after adding isolation plane. Figure 3(c) is the detailed picture. Figure 3(d) is the traditional three-dimensional RG method result. Figure 3(e) is the result after adding isolation plane.

Fig. 3.
figure 3

(a) Geometric diagram of the spatial position of the isolation plane. (b) A slice after adding isolation plane. (c) The detailed picture. (d) The result of traditional method. (e) The result after adding isolation plane.

Full size image

2.2 Confirm the Seed Point

The shape and location of the main pulmonary artery are relatively stable. Therefore, we choose the seed point in the main area of the PA because the area is large and the position is easy to judge. In the region, we select four different seed points, O, P, Q, W (Fig. 3(a)). In order to ensure the reliability of the segmentation, the average of four points is used as the initial value to avoid the inaccuracy of segmentation caused by the contingency of a single point. The method for determining four seed points is detailed in Sect. 2.1.

2.3 Adaptive Three-Dimensional Region Growth

Our improved 3D RG method selects the seeds automatically, and grows according to the criterion after adding an isolation plane. Until there are no pixel points meeting the threshold condition, growth is terminated. Clinically, doctors need to detect subsegmental branches of the PA for PE diagnosis and they hope the branches are conjoint and precise. In order to segment the branches to subsegments, an adaptive parameter is introduced to segment the PA and its branches accurately.

From pulmonary trunk to branches, CT value descends drastically. It is hoped that the finer we segment, the better effect we will get. That is, the seed iterates automatically on the direction in which the CT value is getting lower. Assume that the CT value of initial seed point in pulmonary trunk is \( s_{0} \), and the ending point in the branch is \( s_{n} \). After iterating its 26 neighborhoods in the first loop, the CT value along the direction decreases a bit. For this characteristic, a weight is assigned to the minimum value of 26 points, which makes seed value cuts down gradually. The weight, an adaptive iteration parameter, is called \( \beta \). Through each loop, the CT value of 26 neighborhoods decreases by \( d \), and the entire three-dimensional growth process needs to visit 26 neighborhoods \( n \) times. At the end of the first iteration, seed is alternative to \( s_{1} \) that is used for the next round of growth. Similarly, we obtain \( s_{2} \) to \( s_{n} \). Variables, \( s_{0} ,s_{1} ,s_{2} , \ldots s_{n} \) and \( d \) are known, and calculate \( \beta \). The iteration Eq. (4) is shown below:

$$ s_{n} = s_{0} - \frac{n*d*\beta }{1 + \beta } $$
(4)

Experimental result turns out to be an accurate segmentation by adding an adaptive iteration parameter. The Fig. 4(a) below is a detailed picture of the traditional algorithm. The Fig. 4(b) is a detailed picture of adding adaptive iteration parameter \( \beta \).

Fig. 4.
figure 4

(a) The detailed picture of the traditional algorithm. (b) The detailed picture of adding adaptive iteration parameter \( \beta \).

Full size image

By means of mathematical calculations and experimental verification, the best effect is acquired when \( \beta \) = 0.0015. It is proved that the descent direction of the seed iterated along the desired direction. As a result, the branches of the PA can be segmented to the subsegmental level. The improved method meets the doctor’s requirements because of the accurate branches and vivid reconstruction model. The three-dimensional reconstruction result is in Fig. 5. Figure 5(a) shows the result when \( \beta \) = 0.00015. Figure 5(b) shows the result when \( \beta \) = 0.0015. Figure 5(c) shows the result when \( \beta \) = 0.015. It can be seen that the branch of Fig. 5(b) is finer and more detailed than Fig. 5(a). Figure 5(c) shows that all regions have been grown because the weight to the minimum point is too large. Therefore, when \( \beta \) = 0.0015 (Fig. 5(b)), the result is optimal.

Fig. 5.
figure 5

(a) The result when \( \beta \) = 0.00015. (b) The result when \( \beta \) = 0.0015. (c) The result when \( \beta \) = 0.015.

Full size image

3 Results and Discussion

In the experiment, eight groups of data with similar anatomic characteristics from China-Japan Friendship Hospital are tested. The lung CTPA images are grayscale images. The MPA length of the patients is known, and at the junction of the pulmonary trunk and the heart, the angle between the pulmonary trunk and the XOY plane is about 20°–60°. The CTPA images of the experimental data come from the two devices of the China-Japan Friendship Hospital. The models include GE Healthcare SYSTEMS REVOLUTION CT of General Electric Company and TOSHIBA Aquilion ONE of Toshiba Medical Company. The experiments are carried out using MATLAB2014a on the PC with Inter Core i7-8700 K, CPU 3.70 GHz, and 32 GB RAM. Two results using our method to segment PA of two patients are shown in Fig. 6(a), (b):

Fig. 6.
figure 6

(a) (b) The result using our method of with two patients.

Full size image

Compared with other scholars’ researches, the proposed three-dimensional region growth algorithm has improved in the number of PA branches, branches level, three-dimensional connectivity, and smoothness of pulmonary artery wall. The anatomical characteristic can be simulated veritably to analyze the degree of risk. The paper adopts volume rendering method to draw three-dimensional model. Medically, the level of pulmonary arteries and their branches is defined as follows. For example, the pulmonary trunk is the first level, the LPA is the second level, the left superior lobar artery is the third level, the segment is the fourth level, and the subsegment is the fifth level. Figure 7 shows the results of three other researchers in References [6, 11, 12]. Figure 7(a) shows the lobar arteries and segments and some branches reach the fourth level. But the branches are fractured, and the number of branches is incomplete. The pulmonary veins are mixed with PA, and the artery wall is rough. Figure 7(b) shows the lobar arteries and segments and some branches reach the fourth level. The number of branches is incomplete and the pulmonary artery wall is rough, too. Figure 7(c) shows that the pulmonary vessels are extracted completely, but the PA and pulmonary vein could not be separated.

Fig. 7.
figure 7

(a) (b) (c) Three results of other researchers in References [6, 11, 12].

Full size image

The results of the traditional and improved approach are shown in Fig. 8. Figure 8(a) is the result of traditional method. Figure 8(a) shows the lobar arteries, segments and subsegments. The branches reach the fourth level and the number of branches is relatively complete. The artery wall is slightly rough and the vena cava is mis-segmented. Figure 8(b) is the result of the improved 3D RG approach with adaptive parameter. Figure 8(b) shows the lobar arteries, segments and subsegments. The branches reach the fifth level and the number of branches is complete and accurate. The artery wall is slightly rough and the vena cava is mis-segmented. Figure 8(c) is the result of the improved 3D RG approach adding the isolation plane. Figure 8(c) shows the lobar arteries, segments, and subsegments and the branches reach the fifth level. The number of branches is complete and accurate, and the artery wall is slightly rough. The segmentation effect is confirmed by radiologist that the result is good and achieves the desired goal.

Fig. 8.
figure 8

(a) The result of traditional method. (b) The result of the improved 3D RG approach with adaptive parameter. (c) The result of the improved 3D RG approach adding the isolation plane.

Full size image

4 Conclusions

In this paper, an adaptive three-dimensional region growth method is used to segment pulmonary artery and its branches. On the basis of vector projection theorem and anatomical characteristics, the skeleton and the normal vector are obtained. An isolation plane is introduced to remove the adhesion of the vena cava. The 3D RG algorithm is improved in choosing four seed points automatically and adding an adaptive parameter. Several sets of CTPA data are tested and it is proved that our approach can segment PA and its branches to subsegment level accurately. This method solves the following problems, including the fractured branches, incomplete branches and segmentation leakage. The results meet the requirements of radiologists, moreover, our method lays a solid foundation for the subsequent detection of pulmonary embolism and risk assessment.