Since its introduction in the 1970s, computed tomography (CT) has become an important tool in medical imaging. It is the gold standard in the diagnosis of a large number of different disease entities , and further technological progress has strengthened its diagnostic impact leading to an essential role in clinical practice. To gain full benefit of the increasing resolution of CT images, automatic methods are needed to separate important information from diagnostically irrelevant ones.
Automatic segmentation and analysis of the vessels inside the lung (pulmonary vessels) from CT images is widely used for computer aided diagnosis of vascular diseases , non-rigid image registration , and detection of pulmonary embolism . Our clinical focus is on the detection of pulmonary hypertension (PH), which is a chronic disorder of the pulmonary circulation, marked by an elevated vascular resistance and elevated mean pulmonary artery pressure (mPAP) . Unlike the systematic circulation, the blood pressure in the pulmonary vessels is very difficult to measure. In order to determine pulmonary pressure, an invasive right heart catheterisation must be performed . By the time of diagnosis, PH has usually progressed to late stage and is not reversible any more. By finding a non-invasive way of measuring the pulmonary blood pressure, the number of patients awaiting treatment could be significantly decreased. Therefore, an important aim of clinical PH research is the early diagnosis of pulmonary hypertension.
For determining a measure for pulmonary blood pressure, a segmentation of the blood vessels inside the lung is needed. We present a robust algorithm that combines lung- and airway-segmentation, together with a sophisticated vessel enhancement filter to obtain a proper segmentation of the left and right pulmonary vessel trees, even in patients showing severe pathologies. A schematic of the anatomy of the human lung is shown in Figure 1.
The algorithm is fully automatic, computationally efficient and able to handle large datasets. It was tested on phantom data, the publicly available VESSEL12 challenge dataset, and CT data from 24 patients from a clinical PH study.
1.1 Related Work
Several 3D vessel segmentation algorithms have been presented in the literature up till now. An overview can be found in . Typical approaches include threshold-based algorithms  or fuzzy methods . These approaches have in common that an intensity model is utilized to detect the vessels. Frangi et al. 
presented a technique based on the eigenvalue analysis of the Hessian matrix. This approach was later refined in the popular approaches of and 
, who also take the eigenvectors of the Hessian matrix into account. In
they apply vessel-enhancement, junction-enhancement and nodule-enhancement filters, based on an eigenvector analysis from the gradient vector field. Most of these approaches work well on controls, however, for patients showing pathologies, robust extraction of vascular structures is still an open issue, especially in the case of pulmonary hypertension where vessel pruning occurs.
At the core of our method is a multi-scale vessel enhancement (VE) filter based on the Hessian matrix. It is similar to 
in using the eigenvectors of the Hessian matrix to detect candidate voxels inside the vessels, and computing an offset-medialness boundary measure perpendicular to the estimated vessel direction. The VE response (i.e. medialness) is limited to the right and left lung, which is derived from an intensity-based lung segmentation and a coarse airway tree segmentation with a labelling of left and right main bronchi. After centerline detection from the VE response, a coarse radius estimation of the vessel is obtained using a spherical ray-casting approach. The final segmentation is the output of a globally optimal geodesic active contour model based on total variation and a data term derived from the medialness . Figure 2 shows the flowchart of our automatic vessel segmentation approach.
2.1 Lung and airway segmentation
A prerequisite for our vascular tree extraction is a segmentation of left and right lungs, reactively, to restrict the reconnection of the vessel centerlines. A coarse airway segmentation and labelling of the main bronchi initiates this process. The labeled airway tree is subsequently used to label a coarse, threshold based, lung segmentation. Segmenting the airways also helps in removing false positives of the vessel enhancement filter (see Section 2.2), since the intensity contrast of the airway border and blood vessels is very low, thus leading to incorrect detection of blood vessels at the airway walls.
We automatically detect the airway on the top-most slice of the contrast-enhanced volume, which is a dark circle surrounded by high-intensity tissue, to get a seed point for an iterative 3D region growing algorithm. For the region growing two thresholds are defined: and , where is the CT image, is the seed point and HU denotes a Hounsfield Unit. All connected voxels which fulfill are added to the segmentation. Then the thresholds are updated ( and ) and region growing is restarted with the new segmentation as seed. Convergence of this iterative procedure combines two stopping criteria. We check whether the number of total voxels, or the number of edge voxels of the current segmentation is three times larger than the average of all of the previous voxel counts. In this case a leakage is detected (see Figure 3b), and we produce the final airway tree segmentation with the restricted threshold range from the previous iteration. Figures 3a shows an example of airway segmentation obtained from one patient in the clinical PH study.
After coarse airway segmentation, we perform a left and right lung segmentation to identify a region of interest for later vessel detection. A coarse lung segmentation is obtained using Otsu’s optimal thresholding method . With a connected component analysis, the lung is selected and a 3D hole filling is applied to include vascular structures. The two lungs always merge through the airways, but in some datasets the border between right and left lung is hardly visible, resulting in connected lungs. The airway segmentation is used to separate the coarse lung segmentation. Using a graph representation of the skeleton from the airway segmentation, we detect the carina (where the trachea splits into the left and right main bronchi), and assign different labels to the trachea, right and left bronchi. To label the voxels in the coarse lung segmentation, we calculate shortest paths to the labeled airway tree, thus, splitting it into left and right lung. As a cost function for the shortest path algorithm, we use the gradient magnitude of the CT image and the coarse lung segmentation without the airways, where we give a larger weight to the gradients to prevent connecting via short cuts from touching left and right lungs. We found a fraction of for the gradient weight, and a weight of for the lung segmentation without the airways to perform well in the labeling of our datasets.
As a final step, to remove holes caused by vessels and other high intensity structures inside the lung, the airways are removed from the lung segmentation, and a morphological closing operation is applied several times at each lung separately. We use a six-neighbourhood star-shaped structuring element and 10 closing operations. These parameters remain constant for all datasets. This ensures that the lung segmentations contain the whole lung. The different steps of the lung segmentation can be seen in Figure 4.
2.2 Vessel enhancement
We enhance all vessel-like structures using a modified version of the vessel enhancement (VE) filter proposed by Pock et al. . We extend the VE filter by using multiple radii in addition to a multi-scale pyramid with a down-sampling factor of 1.7. With these parameters, the experiments showed improved discrimination between different vessel radii compared to the widespread use of Gaussian pyramids with a down-sampling factor of two. Further, to be more robust against noise and to get rotation invariant derivatives, we compute our gradients for the boundary measure and for the Hessian matrix using the filter kernels of . The airway- and lung-segmentations from Section 2.1 are used to restrict the vessel enhancement output to the left and right lungs without the airways, respectively.
To get the vessel enhancement filter response, we calculate the eigenvalues , , and (ordered such that ) and the associated eigenvectors , and of the Hessian matrix at each scale . To sort out for bright tubular structures on dark background, we check that and holds. In points that fulfill this condition, the smallest eigenvector gives an estimation for the vessel direction. Perpendicular to the vessel direction, in the cross section plane of the tube given by the eigenvectors and , we evaluate boundary information along circles of different radii . We define the boundary gradient , with
being the CT image convolved with a Gaussian kernel with variance. An initial response is given by the median of the boundary contributions , with , which we denote as . A problem of is, that it also produces responses at isolated edges. To avoid this, a measure of symmetry is introduced:
where is the median absolute deviation of the boundary samples and is the median. The final boundary response is computed as:
To suppress responses at the border of vessels, the gradient magnitude at the center of the vessel is combined with the offset medialness from above:
The final vesselness response
is the maximum response from all different scales and radii . We found 4 scales and radii varying from 1 to 2 pixels, with an increase of 0.3 pixels to have the best performance.
2.3 Centerline extraction
In a non-maximum suppression step inspired by , at each position with a medialness , we sample 8 points on a plane perpendicular to the estimated vessel direction. If the medialness on any of those 8 points is larger than at the current position , the VE response at is set to zero. This leads to disconnected vessel centerline fragments due to branching points, where the tubularity assumption fails. Next, small centerline fragments (less than 5 -connected voxels) are removed, and all maxima lying on the airway border are cleared. To reconnect the centerline fragments, we apply a Dijkstra-like shortest path algorithm. In each lung separately, we connect all centerline candidate points to the center of the heart. As a cost function, we combine the medialness with the gradient magnitudes of the CT image. The separate processing of right and left lung ensures that wrong connections through the mediastinum are avoided. The connected trees of the right and left lung form the final vessel tree.
2.4 Vessel segmentation
At each centerline voxel, 48 points lying on a sphere are sampled from the CT image and the grey values are summed up. This is done for spheres with different diameters. As soon as the normalized sum of grey values drops under 0.6 (empirically found), the radius of the vessel has been found. This gives a more accurate radius estimation compared to the filter radii directly from the VE filter response. This coarse segmentation is then used as an input for a total variation based segmentation , which gives the final segmentation.
3 Experimental results
For testing our algorithm we used the publicly available VESSEL12 challenge dataset as well as 24 datasets from patients who underwent contrast enhanced CT as a part of a clinical PH study at the Ludwig Boltzmann Institute for Lung Vascular Research, Graz. The median size of the CT volumes is pixel. Our hardware consists of an Intel CORE I7-2600K 3.40GHZ with 16 GB RAM and a Geforce 580 GTX with 3 GB RAM.
All of our filters are implemented in CUDA, which is a parallel programming model from NVIDIA, that generates hardware accelerated instructions for NVIDIA graphics processing units (GPUs). Using CUDA we can significantly improve the runtime of our 3D image processing algorithms, due to the high degree of parallelization. Therefore, we can directly work on the full resolution CT data. The runtime of the CUDA implementation of the whole algorithm pipeline ranges from 5 to 10 minutes, whereas a CPU-only implementation needs more than an hour. Limited memory on our GPU restricts the possible size of datasets in the different stages of our pipeline, so we decompose large (i.e. nearly ) CT images into overlapping sub-volumes, which are processed sequentially, with each sub-volume benefiting from the CUDA based parallelization.
3.1 Airway and lung segmentation
On 24 patients we were able to extract the airway tree to at least four generations, which is sufficient for our further processing. Due to GPU based parallelization, one airway segmentation takes on average 69 seconds. Figure 3a show a representative result.
The lung segmentations from the 24 patients included all lung tissue and pulmonary vessels, which has been verified through visual inspection (Figure 5). Processing time for the lung segmentation was on average 132 seconds.
3.2 Phantom data set
We used a phantom of a liver vessel tree, depicted in Figure 6
a, to check the performance of the algorithm and validate its robustness against noise. We successively added Gaussian noise with increasing variance to the phantom data, and calculated the Jaccard index of the ground-truth segmentation with the obtained segmentation. The curve in Figure6b shows how the Jaccard index changes if Gaussian noise with increasing variance is added to the phantom dataset. As long as the variance of the noise is below Hounsfield Units (HU), the performance lies above 93% segmentation overlap.
3.3 VESSEL12 challenge
We applied our algorithm to 20 datasets made available through the VESSEL12 challenge222http://vessel12.grand-challenge.org/. Table 1 shows the results. The algorithm performs very well in terms of specificity, however improvement in the sensitivity is still necessary. This is because the algorithm is optimized for contrast enhanced CT images and for finding vessels even in noisy datasets. Thus vessels smaller than 2mm in diameter are misclassified as noise and are not included in the segmentation. All results from the other participating groups can be found on the official VESSEL12 challenge website.
3.4 Contrast enhanced CT images from our clinical PH study
Our clinical application is the detection of pulmonary hypertension (PH), which is a chronic disorder of the pulmonary circulation, marked by an elevated vascular resistance and elevated mean pulmonary artery pressure (mPAP). Our hypothesis is, that the pulmonary vascular tree shows quantifiable differences between patients with and without PH. One quantifiable property of the vessels is their tortuosity, which is a readout of twistedness . The most common metric of vascular tortuosity is the distance metric (DM), which provides a ratio of the actual vessel length to the Euclidean distance between its endpoints . To determine the tortuosity, the lung vessel centerlines and branching points are extracted. The DM is calculated and compared with the patient’s clinical parameters.
We found a correlation between DM and mPAP of (Spearman correlation coefficient, p<0.01). There was a significant difference between the DM of patients with and without PH (Table 2, p<0.05), thus enabling to discriminate the two groups on our dataset of 24 patients. Two representative vessel segmentation results of the PH study datasets are shown in Figure 7.
We have presented a segmentation approach for vascular structures from contrast-enhanced CT images using a multi-scale vessel enhancement filter and using information from a lung- and airway-segmentation. We achieved very good segmentation results on our 24 patients from a clinical PH study. We also tested the algorithm on non-contrast-enhanced data from the VESSEL12 challenge, where we occupy a midfield position among all participating teams. We see room for improvement in the case of small vessels. Reasons for this performance of the algorithm are the optimization for the contrast-enhanced setup and the use of isotropic CT scans, which is not the case in the VESSEL12 datasets. Due to a parallelized CUDA implementation, our whole vessel tree segmentation and centerline extraction shows a run-time of at most 10 minutes for large CT datasets, without the need for computing on reduced resolutions, thus enabling the potential use in daily clinical routine.
As an important outcome of our work, we showed that tortuosity is correlated with mean pulmonary artery pressure, and our vessel segmentation algorithm can detect the presence of PH. One of the limitations of this study is the small number of patients, which allows only preliminary conclusions. A large scale prospective study to determine the true benefits and constraints of this method is currently in planning. Further, due to the radiation exposure one cannot test the repeatability of the method. This would be necessary to determine its ability for use in disease monitoring and follow-up examinations.
-  G. Agam, S.G Armato, and C. Wu. Vessel tree reconstruction in thoracic ct scans with application to nodule detection. Medical Imaging, IEEE Transactions on, 24(4):486–499, 2005.
-  C. Bauer, T. Pock, E. Sorantin, H. Bischof, and R. Beichel. Segmentation of interwoven 3d tubular tree structures utilizing shape priors and graph cuts. Medical Image Analysis, 14(2):172–184, 2010.
-  E. Bullitt, G. Gerig, S. M. Pizer, W. Lin, and S. R. Aylward. Measuring tortuosity of the intracerebral vasculature from mra images. IEEE Transactions on Medical Imaging, 22:1163–1171, 2003.
-  M. Day and C.J. Topliss. Radiological diagnosis of degenerate change of the subtalar joint: A study comparing the reported degree of osteoarthritis in a plain radiograph when compared with a ct scan. Journal of Bone & Joint Surgery, British Volume, 93-B(SUPP IV):479, 2011.
-  H. Farid and E. P. Simoncelli. Differentiation of discrete multidimensional signals. IEEE Transactions on Image Processing, 13:496–508, 2004.
-  A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever. Multiscale vessel enhancement filtering. In MICCAI, pages 130–137, 1998.
-  N. Galie, M. M. Hoeper, and M. Humbert et al. Guidelines for the diagnosis and treatment of pulmonary hypertension. Eur. Respir. J., 34(6):1219–1263, 2009.
-  J. N. Kaftan, A. Bakai, M. Das, and T. Aach. Locally adaptive fuzzy pulmonary vessel segmentation in contrast enhanced ct data. In ISBI, pages 101–104, 2008.
-  A. P. Kiraly, E. Pichon, D. P. Naidich, and C. L. Novak. Analysis of arterial subtrees affected by pulmonary emboli. In SPIE Medical imaging 2004. Image processing, volume 5370, pages 1720–1729, 2004.
-  K. Krissian, G. Malandain, N. Ayache, R. Vaillant, and Y. Trousset. Model-based detection of tubular structures in 3d images. Computer Vision and Image Understanding, 80(2):130–171, 2000.
-  T. J. Kulik, J. E. Harris, and D. B. McElhinney. The impact of pulmonary venous hypertension on the pulmonary circulation in the young. Congenit Heart Dis, 6(6):603–607, 2011.
-  D. Lesage, E. D. Angelini, I. Bloch, and G. Funka-Lea. A review of 3d vessel lumen segmentation techniques: Models, features and extraction schemes. Medical Image Analysis, 13(6):819–845, 2009.
-  S. Maeda, H. Kim, J. K. Tan, S. Ishikawa, S. Murakami, and T. Aoki. Nonrigid image registration method for thoracic ct images using vessel structure information. In SMC, pages 1413–1417, 2012.
-  N. Otsu. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man and Cybernetics, 9(1):62–66, 1979.
-  T. Pock, R. Beichel, and H. Bischof. A novel robust tube detection filter for 3d centerline extraction. In Image Analysis, volume 3540 of Lecture Notes in Computer Science, pages 481–490. Springer Berlin Heidelberg, 2005.
C. Reinbacher, T. Pock, C. Bauer, and H. Bischof.
Variational segmentation of elongated volumetric structures.
IEEE Conference on Computer Vision and Pattern Recognition, CVPR, pages 3177–3184, 2010.
-  H. Shikata, G. McLennan, E. A. Hoffman, and M. Sonka. Segmentation of pulmonary vascular trees from thoracic 3d ct images. Journal of Biomedical Imaging, 24:1–11, 2009.
-  C. Zhou, H. P. Chan, B. Sahiner, L. M. Hadjiiski, A. Chughtai, S. Patel, J. Wei, J. Ge, P. N. Cascade, and E. A. Kazerooni. Automatic multiscale enhancement and segmentation of pulmonary vessels in CT pulmonary angiography images for CAD applications. Med Phys, 34(12):4567–4577, 2007.