Automatic segmentation of arterial tree from 3D computed tomographic pulmonary angiography (CTPA) scans

Abstract Pulmonary embolism (PE) and other pulmonary vascular diseases, have been found associated with the changes in arterial morphology. To detect arterial changes, we propose a novel, fully automatic method that can extract pulmonary arterial tree in computed tomographic pulmonary angiography (CTPA) images. The approach is based on the fuzzy connectedness framework, combined with 3D vessel enhancement and Harris Corner detection to achieve accurate segmentation. The effectiveness and robustness of the method is validated in clinical datasets consisting of 10 CT angiography scans (6 without PE and 4 with PE). The performance of our method is compared with manual classification and machine learning method based on random forest. Our method achieves a mean accuracy of 92% when compared to manual reference, which is higher than the 89% accuracy achieved by machine learning. This performance of the segmentation for pulmonary arteries may provide a basis for the CAD application of PE.


Introduction
Pulmonary vascular diseases often occur in pulmonary arteries, such as pulmonary emboli (PE), pulmonary arterial hypertension (PAH) or pulmonary aneurysm, featured by the changes of arterial morphology, which can be captured by computed tomography (CT) scanners. Recently, CT pulmonary angiography (CTPA) has been regarded as the gold standard for the diagnosis of PE in clinical practice, as its advantages over both traditional angiography and ventilation perfusion scans.
But the large amount of CT images and imaging noise make it time-consuming and difficult for physicians to diagnose PE from CTPA images. Therefore, the automatic segmentation and visualization of pulmonary arteries is a critical step for the disease diagnosis [1]. Because of the multiple connections between arteries and other anatomical structures, such as the pulmonary veins intertwining with the pulmonary artery, segmentation of the pulmonary arterial tree is a complex problem. Even in contrastenhanced scans, sometimes it is still difficult to distinguish arteries and veins robustly, if the imaging delay from the contrast injection is not controlled properly [2].
Although a large amount of work has been devoted to the segmentation of vessels in different parts of body, the automatic segmentation of pulmonary arteries has not been performed successfully, in respect of the segmentation performance and automaticity. The approaches of pulmonary artery segmentation are normally built on the basis of vessel enhancement or vessel segmentation. The performance is limited because these methods are often interfered by the pulmonary vein. Most 3D pulmonary artery segmentation methods take structural and anatomical information into account, either tracking vessels starting at given seed points or calculating a voxel-wise distinction of arteries and veins [3][4][5], and incorporating different anatomical features like proximity of arteries and bronchi [6]. Saha et al. [7] introduced fuzzy distance transform to realize the separation of arteries and veins, however, in which the interactive refinement is complex and indispensable. Christian et al. [2] used fully automatic method based on integer program to separate arteries and veins in chest CT images, but it may cause the discontinuity of arteries which have an influence on the diagnosis of pulmonary embolism.
To resolve these problems, we introduce a novel, fully automatic approach which segment arteries from volumetric chest CTPA images. We were inspired by the anatomic structure of the pulmonary artery, automatically tracked the directions of artery which starts from pulmonary trunk. There are no interactive options in the method, the degree of automation is comparative to machine learning [8,9].

Method
Our work for pulmonary artery segmentation is shown in Figure 1. We started with lung segmentation based on 3D region growing from CTPA images. Next, the 3D multi-scale vessel enhancement filter based on Hessian matrix was applied to enhance vessels. For later pulmonary artery segmentation, the Harris Corner detection was applied to remove the vena cava, which are adjacent to pulmonary artery trunk in CT images. Finally, at the core of our approach is a fuzzy connectedness framework, by computing the probability that a voxel belongs to the artery tree.

Lung segmentation
In order to improve the performance of artery segmentation, the lung segmentation is a prerequisite for our algorithm. The pulmonary region was automatically extracted by 3D region growing algorithm. For the region growing, a threshold is defined as: where I is the CT image, x 1 and x 2 are two points with maximum and minimum intensity connected to the seed point on the cross section. All 26-neighbourhood connected voxels which fulfill jI(x) ÀI seed j < th are added to the segmentation until no new points added. Then, to remove holes caused by vessels and other high intensity structures inside the lung, 3D hole filling and morphological close operation are applied [10].

Vessel enhancement
To eliminate the effects of lymphoid tissues, the vessel enhancement filter based on the 3D Hessian matrix was used to differentiate vascular tree from other structures, according to the second order gradient information of images. In addition, multi-scale Gaussian filter was used to discriminate vessels with different radius. In order to get the vessel enhancement filter response at voxelr ¼ x; y; z ð Þ in a 3D image, we calculated the eigenvalues k 1 , k 2 and k 3 (ordered such that jk 1 j > jk 2 j > jk 3 j) and the associated eigenvectors are e 1 , e 2 and e 3 of the Hessian matrix at each scale r s [10]. The multi-scale response function: where n is a constant, plays an important role in enhancement of the tubular and the blob-like structures in Equation (2). For tubular structure of vessels branch, the eigenvalues subjected to jk 1 j % jk 2 j, jk 1 j ) jk 3 j, and jk 3 j % 0. And the vessel bifurcation which forms a blob-like structure when the vessel splits into two or more branches, with the eigenvalues jk 1 j % jk 2 j % jk 3 j. While the lymphoid tissues surrounding the pulmonary vessels which generally are plate-like structures with jk 1 j ) jk 2 j, and jk 2 j % jk 3 j % 0 [11]. The response functions of blob-like and tubular structures in images were enhanced, and plate-like structures were suppressed, achieving aims of vessels enhancement.

Harris corner detection
Due to the contrast injection from vena cava, the vena cave greatly interfered the segmentation of pulmonary artery, because the vena cava has high image intensity and is close to the pulmonary artery trunk. To solve this problem, Harris corner detector was used to detect the boundaries between vena cava and pulmonary artery trunk, which depends on the autocorrelation of gradient magnitude image to analyze the corner points in image [12]. The pixels have been classified on the basis of the relation of eigenvalues as follows, both high eigenvalues representing the corner pixel, while both small eigenvalues indicating that pixel is in flat region, and only one eigenvalue representing the edge pixel. The corner points between vena cava and pulmonary artery trunk were detected by setting a threshold for Harris detector, and then fitted to a line to separate vena cava from pulmonary artery trunk through 3D region growing method.

Fuzzy connectedness
The 3D fuzzy connectedness algorithm was introduced to the automatic segmentation of pulmonary arterial tree, which is based on a global fuzzy relation that assigns a strength of connectedness to every pair of voxels in an image [13]. Given an object O with a seed point s and its background B are classified by dividing the set of voxels in the volumetric image in a way of computing its affinity m k with the seed point. The "strength of connectedness" of two distant voxels c and d along a certain path p(c, d) within the image is simply the smallest pairwise local affinity along this path. If two pixels c and d are adjacent, their affinity m k (c, d) decides whether they belong to the same class or not. The adjacency component ma is a non-increasing function of the distance in voxels jjc À djj. One general form is where f(c) and f(d) correspond to the CT intensities at pixels c and d respectively, h 1 and h 2 are the Gaussian metric which is a function of mean and standard deviation of the CT intensity values of the object, and x 1 þ x 2 ¼ 1 [14].
If the affinity is larger than a certain threshold, the respond voxels will be defined as object. Through several iterations, the strength of relative paths can be weaken, the responsive voxels are excluded from the class of object. The advantage of our method is the automaticity that even the only one seed point is chosen automatically at the point of highest intensity, owing to our images are contrast enhanced on the pulmonary artery.

Data preparation
For testing our algorithm, we used the clinical datasets which obtained from China-Japan friendship hospital. The thoracic regions were imaged using Toshiba Medical Systems multi-detector CT scanner at 120kVp under automatic exposure control. Contrast materials were injected into patients at 3.0 ml/sec for obtaining contrast data. The dataset consists of 6 non-embolism and 4 embolism patients. All CTPA images were read and labelled by trained radiologists through 3D Slicer.

Experimental design
In section of lung segmentation, when consider the connectivity between the lung and the airway, an initial seed point was set automatically inside the airway, which had low CT value and near the middle of the image of the upper chest. The results for lung segmentation were shown in Figure 2. In vessel enhancement filter, n ¼ 0.7 may achieve the differences of tubular-like and blob-like structures from plate-like structures [11]. Figure 3 showed the effect of the parameter n on the vessel enhancement, indicating that n was not so sensitive, since the results of vessel enhancement did not vary significantly when n ranged from 0.5 to 0.7. The scales of Gaussian filter was set from 1 to 12 correspond to different vessel sizes according to the image size. The results for vessel enhancement was shown in Figure 4, where the pseudo-color showed the likeliness of enhanced vessels, the red parts of which were obviously enhanced. The Harris Corner detection was introduced, and the result for corner detection was shown in Figure 5. The cross section between pulmonary artery and vena cava was highlighted in different colors.
In fuzzy connectedness algorithm, we set the threshold of m k (c, d) as 0.1, which distinguished the voxel with the fuzzy connectedness smaller than 0.1   as the same object as the seed point [13]. The influence of the threshold on the results of pulmonary artery segmentation was shown in Figure 6. And h 1 and h 2 are automatically computed by 3 Â 3 Â 3 volumes round the voxel. The final results were shown in Figure 7, where the 3D surface rendering was used to reconstruct the pulmonary artery tree. The patients (g)-(j) suffered with PE, and they are different from (a)-(f) in vessel density and distribution, indicating that our method could be helpful to distinguish PE. The average processing time of each case is about 3.5 minute on a machine with Intel Core (i5) CPU, 8 GB RAM and Windows 7 OS.

Quantitative evaluation
In this part, accuracy of our method is defined by comparing its results Artery seg with manual labeling of arteries in Artery label , which was regarded as ground truth. Here, we have defined segmentation performance measures in terms of voxel-based agreement (Dice coefficient) Agreement voxel , the number of vessels on each slice Similarity number , Sensitivity voxel and Precision voxel as follows: Artery seg þ Artery label j j (4) Artery label (6) Precision voxel ¼ Artery seg \ Artery label Artery seg (7) where Agreement voxel was used to measure reproducibility between automatic and ground truth  segmentations. Sensitivity voxel denotes the percentages of the volume that vessel in the ground truth was correctly classified in the manual segmentation, and Precision voxel was used to calculate the percentage of the vessel in the ground truth that was correctly classified in the segmentation results. Table 1 shows the results to evaluate our method and machine learning method based on random forest proposed by Zhao [8]. The machine learning method consists of 3 steps. Firstly, multi-scale representations of images are obtained via the Gaussian pyramid. Then, a large number of patches are randomly selected from multiscale representations and a sparse auto-encoder is trained using these patches. Finally, the random forest method is applied to segment lung arteries with features and labelled data. The set of the parameters in this method is as follow: (1) the level of Gaussian pyramid is 6; the window size is 5 Â 5; (2) the patch size is 5 Â 5; the layer of sparse auto-encoder is 3; the number of neurons in the input and output layers is 25; the number of neurons in the hidden layer is 32; (3) the number of trees in random forest is 500 and the number of randomly selected feature is 8. From Figure 8, we can find that the pulmonary arteries segmented by our method may have less area than the ground truth, namely the manual segmentation results, in some cases, which may be caused by the border-effect in the vessel enhancement. The segmentation results from previous conventional randomforest-based method random forest algorithm showed some discontinuity of vessels marked in Figure 3, which was also reported by former researchers in the studies of machine learning algorithm. To summarize the results, the algorithm performs very well in voxelbased agreement, vessel numbers agreement and precision, however improvement in the sensitivity is still necessary. From Table 1, we can see that the proposed method achieves a mean voxel-based agreement with manual references of 92.16%, and the precision is 94.32%, which are both significantly higher  than machine learning. Our method achieves a mean number similarity of 90.97%, which is dramatically greater than machine learning of 80.73%. On the other hand, the mean value of sensitivity of our method is slightly lower than machine learning. To conclude, for majority of datasets we get an excellent performance than the random-forest-based method.

Discussion and conclusion
Our method was based on the anatomical structure of lung, by tracking the path of arteries or contrast agent inside lung to segment arteries. The novelty of our work is that we proposed an improved fuzzy connectedness algorithm combined with Harris corner detection and Hessian filter, which achieved fully automatic arterial tree segmentation, and could be used in the CAD application of PE. Hessian filter and Harris corner detection was used to separate arteries from veins and other organs around heart. In our method, there are few parameters need to be adjusted, which is easier to form into an automatic image-based CAD system, than other methods like level set. For most cases, the automaticity of our method is comparable to the machine learning method.
Although a few promising methods have been proposed for arterial segmentation, they made strict demands on the imaging process, which made it difficult to compare these results with ours [15,16]. When comparing our method to machine learning, the automaticity of the two methods is comparable. More importantly, the results and connectivity of our method is better in small datasets, although the performance of machine learning may improve with the increasing of datasets. Therefore, more subjects are needed to add to the datasets to validate the effect of our method. For now, our datasets are proved to be challenging due to the complex clinical samples. In the patients with PE, the structure of pulmonary arteries may be irregular, less like a tubular structure, which may reduce the accuracy of the segmentation. However, for the images from Patient (g) -(i), who suffered from PE, the performance of our method is as good as that from the normal. For Patient j, the segmentation results are relatively low, probably due to the presence of bronchiectasis. In the future, more pulmonary embolism datasets should be tested through our method. Since PE reduces the arterial distribution in CT images, our method is expected to calculate the artery density or distribution, which could be used for the aided detection of PE.

Disclosure statement
No potential conflict of interest was reported by the authors.  . The comparison between the segmentation of pulmonary arteries by the ground truth, our method and previous conventional random-forest-based method. The images are from subject a. The first column is the original image. The second is the manual segmentation results of pulmonary arteries, while the third and fourth column is the results of our method and the random-forest-based method respectively.