Assessment of fully automatic segmentation of pulmonary artery and aorta on noncontrast CT with optimal surface graph cuts.

PURPOSE: Accurate segmentation of the pulmonary arteries and aorta is important due to the association of the diameter and the shape of these vessels with several cardiovascular diseases and with the risk of exacerbations and death in patients with chronic obstructive pulmonary disease. We propose a fully automatic method based on an optimal surface graph-cut algorithm to quantify the full 3D shape and the diameters of the pulmonary arteries and aorta in noncontrast computed tomography (CT) scans.
METHODS: The proposed algorithm first extracts seed points in the right and left pulmonary arteries, the pulmonary trunk, and the ascending and descending aorta by using multi-atlas registration. Subsequently, the centerlines of the pulmonary arteries and aorta are extracted by a minimum cost path tracking between the extracted seed points, with a cost based on a combination of lumen intensity similarity and multiscale medialness in three planes. The centerlines are refined by applying the path tracking algorithm to curved multiplanar reformatted scans and are then smoothed and dilated nonuniformly according to the extracted local vessel radius from the medialness filter. The resulting coarse estimates of the vessels are used as initialization for a graph-cut segmentation. Once the vessels are segmented, the diameters of the pulmonary artery (PA) and the ascending aorta (AA) and the P A : A A ratio are automatically calculated both in a single axial slice and in a 10 mm volume around the automatically extracted PA bifurcation level. The method is evaluated on noncontrast CT scans from the Danish Lung Cancer Screening Trial (DLCST). Segmentation accuracy is determined by comparing with manual annotations on 25 CT scans. Intraclass correlation (ICC) between manual and automatic diameters, both measured in axial slices at the PA bifurcation level, is computed on an additional 200 CT scans. Repeatability of the automated 3D volumetric diameter and P A : A A ratio calculations (perpendicular to the vessel axis) are evaluated on 118 scan-rescan pairs with an average in-between time of 3 months.
RESULTS: We obtained a Dice segmentation overlap of 0.94 ± 0.02 for pulmonary arteries and 0.96 ± 0.01 for the aorta, with a mean surface distance of 0.62 ± 0.33 mm and 0.43 ± 0.07 mm, respectively. ICC between manual and automatic in-slice diameter measures was 0.92 for PA, 0.97 for AA, and 0.90 for the P A : A A ratio, and for automatic diameters in 3D volumes around the PA bifurcation level between scan and rescan was 0.89, 0.95, and 0.86, respectively.
CONCLUSION: The proposed automatic segmentation method can reliably extract diameters of the large arteries in non-ECG-gated noncontrast CT scans such as are acquired in lung cancer screening.

INTRODUCTION

Cardiovascular diseases and chronic obstructive pulmonary disease (COPD) are among the major leading causes of death globally. In the search for early identification of individuals at risk of cardiovascular disease in COPD, , imaging‐based assessments of the shape and size of the pulmonary artery (PA) and aorta have rapidly gained interest. Changes in these two large arteries may indicate cardiovascular diseases, including pulmonary hypertension, , aortic dilatation and aortic aneurysm, and coarctation of the aorta. The PA to ascending aorta (AA) diameter ratio () at the level of PA bifurcation is shown to be associated with the presence of pulmonary arterial hypertension and is associated with poorer health status increased risk of severe exacerbations, , and increased mortality in patients with COPD.

Performing diameter measurements manually is labor‐intensive and time‐consuming, and has high intra and interobserver variability. Diameter measurements derived from 3D segmentations are more reliable but are even more time‐consuming to obtain manually. To accurately assess the pulmonary arteries and aorta, automatic 3D segmentation is, therefore, desirable.

With the growing use of low‐dose noncontrast thoracic computed tomography (CT) scans for lung cancer screening, , , there is an opportunity to measure the pulmonary arteries and aorta in these scans in order to investigate the presence of early‐stage cardiovascular disease and/or predict complications in patients with COPD. However, in noncontrast CT, segmentation of the aorta and especially the PA is challenging due to their proximity to other structures with similar intensity values. In nonelectrocardiography (ECG)‐gated CT as is commonly used in lung screening, additional challenges are motion artifacts and unclear vessel boundaries at the regions close to the heart.

In the literature, automated segmentation methods of the pulmonary arteries and aorta have been presented mainly for magnetic resonance imaging and contrast‐enhanced CT angiography (CTA), , , which have high contrast between vessels, fat, and surrounding muscles. However, these methods do not translate well to noncontrast CT, where the vessel boundaries are not well defined in many places. Even though MRI and CTA have high contrast, CT is widely used in clinical practice. Compared to MRI, CT has a faster acquisition time, better isotropic spatial resolution, convenience, and easier access. Furthermore, the iodinated contrast material in CTA can produce undesired side effects, such as allergic reactions and kidney damage. Therefore, it is desired to avoid contrast when the vasculature is likely to be visible on noncontrast‐enhanced CT. Many patients with COPD or at risk of developing cardiovascular disease undergo a low‐dose, non‐ECG‐gated, noncontrast thoracic CT for lung cancer screening. The ECG gating process requires highly overlapping slices, which exposes patients to a higher radiation dose. Therefore, to reduce the risk of radiation exposure, it is common to perform low‐dose, non‐ECG‐gated, noncontrast CT for lung cancer screening.

Despite the widespread use of noncontrast CT in clinical practice, relatively fewer studies can be found on these scans, on the segmentation of the aorta , , , , , , , , , , , , and especially the PA containing the pulmonary trunk, left and right pulmonary arteries. , ,

Among existing segmentation methods, those using a shape prior , , , generally obtain good segmentation results on noncontrast CT scans. Xie et al. , employed a cylinder matching method to extract the centerline of the PA trunk and aorta. To segment the vessels, they used geometric constraints from adjacent organs obtained from a precomputed anatomy label map. Although the obtained results are good, an anatomy map is not always available.

Another approach is to use graph‐cut methods with shape priors. Graph cuts can achieve a global optimum with low processing times and it is possible to incorporate shape constraints in the graph structure. In the approach presented in Ref. , Boykov et al. showed that a convex cost function of a combinatorial problem described as a graph, an optimal solution is obtained by a graph‐cut. As an optimal surface graph is a type of Boykov style graph, then the condition holds, and as the graph is constructed to be convex, then optimality should be guaranteed. Graph‐cut methods have been applied to different imaging modalities for artery segmentation , and have obtained promising results in many tasks. Deng et al. proposed a graph‐cut method using random forest‐based discriminative features on noncontrast CT for aorta segmentation. They achieved a high segmentation performance in the abdominal aorta; however, they have not applied their method to the PA, which is more challenging to segment.

Deep learning‐based techniques are being used widely recently and have shown satisfying results. Noothout et al. proposed a dilated convolutional neural network (CNN) for aorta segmentation that classifies the voxels into axial, coronal, and sagittal CT slices and subsequently obtains the final segmentation result by averaging the resulting probabilities over the three directions. Haq et al. trained a multilabeled network based on the DeepLab architecture and segmented both the PA and aorta on noncontrast CT scans. For improved segmentation results on noncontrast CT, Bruns et al. and Lartaud et al. trained their multilabel (heart and large arteries) deep learning‐based method with virtual noncontrast CT scans generated from dual‐energy CTA scans. Bruns et al. segmented the AA and pulmonary trunk by using a CNN architecture, and Lartaud et al. segmented the thoracic aorta by using a six‐label U‐Net. Despite the high per‐voxel accuracy of the deep learning‐based methods, preserving the vessel topology is a challenge in these methods. Therefore, deep learning‐based methods are not guaranteed to result in smooth structures with vessel topology.

In this paper, we present an optimal surface graph cut‐based method that enforces topology constraints and incorporates a shape prior in the graph structure to segment the pulmonary arteries and aorta. We adopt the optimal surface graph‐cut approach by Petersen et al., originally proposed for airway segmentation, which incorporates a shape prior via constructing the graph based on flow lines traced from an initial, smoothed segmentation. The nonintersecting flow lines guarantee nonself‐intersecting surfaces and make it possible to segment high curvature areas, such as the bifurcation of the PA and the aortic arch while guaranteeing a shape that is similar to the initialization shape.

A preliminary version of this work is presented in Ref. . In the current work, the proposed method is fully automated and includes a multi‐atlas registration technique to automatically extract seed points and a landmark for the level of the PA bifurcation. At the level of the extracted landmark, the PA and AA diameters and ratio are automatically extracted. The main contribution of the current work is in the extensive validation of the accuracy, robustness, and reproducibility of the optimal surface graph‐cut method in the challenging segmentation problem of the PA and aorta segmentation in noncontrast CT using a large dataset. An additional contribution is in fully automating the process of seed point, centerline, diameter measurement, and ratio extraction. Unlike our previous work focusing on aorta segmentation, in this paper, we segment both the PA and aorta with a more robust landmark detection technique.

METHODS

The main steps of our proposed method are (1) preprocessing; (2) automatic seed point and landmark extraction with multi‐atlas registration; (3) centerline extraction for vessel localization; (4) vessel segmentation using an optimal surface graphcut algorithm; (5) 3D diameter measurement and biomarker extraction. An overview of our method is shown in Figure 1.

FIGURE 1

An overview of the method. Abbreviations: AA, ascending aorta; DA, descending aorta; LPA, left pulmonary artery; PABIFL, pulmonary artery bifurcation level; PAT, pulmonary artery trunk; RPA, right pulmonary artery

Preprocessing

To reduce the unnecessary computational cost, a bounding box is calculated around the lungs. The lungs are segmented using thresholding and morphological smoothing similar to Lo et al. Thereafter, the scans are cropped using a bounding box around the lungs Figure 2a. Then, low‐ and high‐intensity values are clipped, to prevent the centerline extraction and vessel segmentation attracting to the heart–lung or bone–lung borders, which often have a higher gradient than the boundary of the vessels of interest. Truncating intensities higher than 150 HU, such as presented in bones, or lower than –150 HU, generally presented in lungs, makes the vessel borders relatively stronger and easier to detect. The gradient magnitude of an axial slice shown in Figure 2b–e illustrates the effect of preprocessing on enhancing the edges of the pulmonary arteries and aorta.

FIGURE 2

An axial view of a CT scan with the bounding box in red is in (a), zoomed in view of the original scan in (b), and the preprocessed scan in (d). The corresponding gradient magnitude of the original scan is in (c) and the gradient magnitude of the pre processed scan is in (e). Strong edges at lung and bone borders (b, c) are removed during preprocessing (d, e). In the second row (f–i), images are overlaid with the manual PA (green) and aorta (yellow) annotations

Seed point extraction

Seed points are obtained with a multi‐atlas registration method similar to the one presented in Tang et al. A set of atlas images and their corresponding label images , are used. includes manual seed points at the left and right pulmonary arteries before the secondary bifurcations, PA trunk, PA bifurcation point, AA at the sinotubular junction, and descending aorta at the diaphragm level. To address the large variation in the shape and size of these arteries, is set to 25 atlases.

The multi‐atlas registration approach consists of three stages. First, each of the atlas images is registered to the target image using an affine transformation followed by a nonrigid registration using a B‐spline transformation model. Normalized mutual information is used as the similarity metric. The 10 registered atlas images with the highest similarity (also defined by normalized mutual information) to are selected, and the corresponding transformations are applied to images to propagate the seed points to . Finally, deformed images are combined and the final seed points are obtained by averaging the seed point locations per label.

The seed points of the aorta, left, right, and the trunk of the PA are used in the next step to initialize the centerline extraction. The seed point of the PA bifurcation level is the landmark level for measuring the PA and aorta diameters as well as the ratio.

Centerline extraction

We applied a minimum path tracking algorithm to extract the vessel centerlines between the automatically extracted seed points. For the cost function, we used a weighted combination of a medialness filter, , and a lumen intensity similarity filter , similar to Tang et al. Medialness, , is a multiscale filter which uses the circularity assumption and accumulates edge responses along different circle sizes, within a defined radius range of . This gives strong responses in the voxels in the center of the vessel and drops rapidly toward the vessel boundary. Simultaneously, we extracted an estimate of the vessel radius at each voxel, by extracting the radius of the circle with the strongest edge response at each voxel location. We applied the medialness filter in the axial, coronal, and sagittal planes, with a different radius range for the pulmonary arteries and the aorta. The medialness for the pulmonary arteries is defined to be the maximum medialness of axial, coronal, and sagittal medialness filters and for the aorta, it is defined to be the maximum medialness of axial and coronal medialness filters. Due to the anatomical orientation of the aorta, the descending and AA are mainly circular in the axial plane, and the aortic arch is mainly circular in the coronal plane. Therefore, only the axial and coronal medialness filters are used for the aorta, and the sagittal plane is discarded from the computation.

In regions with unclear vessel boundaries, medialness alone is not sufficient to ensure correct centerlines. To prevent the vessel centerline from moving outside the vessel lumen, we added a lumen intensity similarity term defined as: where is the intensity of the voxel at position . Tang et al. defined and as the mean and standard deviation of intensity in small regions around the seed points. In our case, they are selected as the mean and standard deviation of intensity in the manual annotations of the pulmonary arteries and aorta in 25 CT scans, which is HU and HU. The lumen intensity similarity makes the structures with higher intensity, such as bones, or lower intensity, such as lungs, get a low response. Although selecting a fixed value for and based on the intensities of the manual annotations may make the method more dependent on the training scan protocol, it ensures that the centerline will stay in the vessel even if the seed point is extracted close to the boundary of the vessel. In cases where the seed points are extracted close to the vessel boundary, calculating and in a local region around the seed point as in Tang et al. can result in wrong estimates and might misdirect the centerline to adjacent structures with slightly different intensities.

From both and , high responses were obtained in the vessel center and low responses in the background. The cost function was defined by an inverted combination of weighted and , where the factors and control the importance of each term, respectively, (). From the constructed cost function, the minimum cost path was obtained by applying Dijkstra's algorithm between the automatically extracted seed points for each vessel, that is, one path between the endpoints of the aorta and two paths between the two endpoints of the left and right pulmonary arteries and its trunk.

Finally, to improve centerline accuracy in areas with high curvature, the centerlines were smoothed with a Gaussian filter with standard deviation and then were refined by recomputing the minimum cost path after curved multiplanar reformatting (CMPR) perpendicular to the previous centerline. An improved estimate of vessel radius was extracted from medialness at CMPR step. The centerline for the entire PA was obtained as the union of centerlines for the right and left pulmonary arteries.

Vessel segmentation

To segment the vessels, we applied an optimal surface graph‐cut algorithm that finds a globally optimal solution of a given cost function. By using nonintersecting columns based on flow lines from a predefined initial shape, (self)intersecting surfaces can be avoided and the topology of the prior shape is preserved. This makes segmentation of high curvature surfaces possible.

To construct a graph with nodes and edges, an initial coarse segmentation is used to generate the graph columns. We used a nonuniform morphological dilation of the vessel centerlines as the initial segmentation. Each centerline point was dilated with a spherical structuring element, with a radius extracted from the radius map . This nonuniform centerline dilation provides information about the shape of the vessel and results in a more accurate vessel surface than the uniform centerline dilation used in Arias et al. In the rest of the paper, we will refer to this dilatation as the coarse initial segmentation. Once we computed the coarse initial segmentation, we converted it to a mesh and constructed graph columns based on flow lines. To ensure sufficient boundary accuracy and to help make the graph smoothness penalties work evenly across the surface, the resolution of the initial mesh was set isotropically and at higher resolution than the original image to . The flow lines were obtained by tracing the gradient of a smoothed version of the initial segmentation. A Gaussian kernel with a standard deviation was used to smooth the initial segmentation. The flow lines were traced from each mesh node, inward in the gradient direction and outward in the negative gradient direction. Graph nodes were sampled in graph columns, at regular arc length intervals along the flow lines with a sampling interval set to . The graph nodes represent the possible image positions the vessel surface can take. A set of edges, which represent the association between nodes, connects the nodes in the graph. The two consecutive nodes in the same column are connected by directed edges named intracolumn edges. The cost of these edges represents local image information associated with the border location and is chosen as the first‐order derivative of the image intensity along the graph column. In CT scans, the vessels have a higher image intensity than the background, therefore, the intensity transitions from high to low and the cost gets to its minimum in the border.

To encourage a smooth segmentation similar in shape to the initialization, “smoothness penalty” edges (penalizing nonsmoothness) were added as in Ref. . These intercolumn edges connect the nodes in adjacent columns with a constant cost , penalizing solutions that deviate from the original shape. A minimum graph cut, minimizing the total cost of edges being cut and separating the graph vertices into vessel and background, provided the final segmentation. This minimization was solved with a min‐cut/max‐flow optimization algorithm. Details of the optimal surface graph segmentation approach can be found in Petersen et al.

Diameter measurement and PA:AA ratio

In clinical practice, the diameter of the PA and AA, and the ratio are measured manually, usually in an axial slice in the level of PA bifurcation. This level is defined as the axial slice where the pulmonary trunk bifurcates, ideally where the right and left pulmonary arteries appear to be of similar size. The selected slice may be one of a few axial slices that fit this criterion, which is likely to lead to inconsistent measurements between annotators and also between baseline and follow‐up scans. Therefore, to ensure a consistent measurement, we calculated the diameters from a 3D vessel segment of 10 mm length from the segmented PA and aorta. The segment is selected perpendicular to the vessel centerline from 5 mm before to 5 mm after the automatically extracted landmark for the level of the PA bifurcation. Assuming a circular cross‐section and a negligible curvature, the average diameter is computed from this segment as , where the 0.1 factor is applied to average the diameter of the cylindrical shape with 10 mm height. Subsequently, from the average diameters, the ratio is calculated. A 3D view of the segmented pulmonary arteries and aorta with the 10 mm segment and the cross‐sections in the extracted landmark level is shown in Figure 1.

For a meaningful comparison with diameters that were measured manually in axial slices, the automatic diameters were also calculated in the axial slice at the same level and in the same direction as the manual diameter measurement. Subsequently, the automatically extracted ratio in the axial view was compared to the manually measured ratio in the same view.

EXPERIMENTS

Dataset and manual annotation

The 471 CT scans used in this study are from the Danish Lung Cancer Screening Trial (DLCST). The study was approved by the Ethical Committee of Copenhagen County and funded by the Danish Ministry of Interior and Health. A multidetector CT scanner ( IDT 16 row scanner, Philips Medical Systems) was used to acquire scans at at maximum inspiration breath‐hold and without cardiac gating. The scans have an in‐plane isotropic resolution of and 1 mm slice thickness. Participants were current or former smokers between 50 and 70 years of age. Baseline CT scans of 235 participants, randomly selected, were used for parameter optimization, method development, and method evaluation. In addition, to assess the repeatability of the proposed method, 118 additional participants were selected who had a baseline scan and a repeat scan after on average 3 months [minimum 2, maximum 5]. Table 1 describes the data used for each experiment.

TABLE 1

Number of CT scans used for each experiment

Experiments	CT scans	CT in atlas	Manual annotation
Centerline parameter optimization	10	–	Centerline and seed points
Method validation with full 3D annotations	25	24	3D volume and centerline
Qualitative assessment with visual inspection	436	25	–
In‐slice diameter and PA:AA ratio assessment	200 out of 436	25	2D in‐slice diameter
Repeatability and 3D diameter assessment	236 out of 436 (118 scan–118 rescan)	25	–

Note: In total, 471 CT scans are used for this study.

Manual annotations of the centerline and contours of the pulmonary arteries and aorta were made using an in‐house annotation tool developed in MeVisLab* with a similar framework as was described previously for carotid artery segmentation.

With this tool, first, the centerline of the vessels was drawn manually using the axial, coronal, and sagittal views, to compute a CMPR. Then, the centerlines were checked and modified if needed in the axial view of the CMPR generated every 1 mm along the centerline. In longitudinal views at six different angles, equally spaced every , longitudinal contours of the vessel were drawn manually. Subsequently, cross‐sectional contours were computed using spline interpolation through the intersection points of the longitudinal contours with the cross‐sectional planes. Finally, after checking the cross‐sectional contours in all cross‐sections and adjusting them if required, the contours were converted to a 3D binary image using variational interpolation. For all cases, the window level/width for annotation was set to 200 HU/600 HU.

For the pulmonary arteries, a centerline and binary segmentation image were first created for the pulmonary trunk + left PA and trunk + right PA individually. The segmentation for the entire pulmonary was obtained as the union of these two segmentations.

With this tool, centerlines of 35 CT scans (10 for optimizing the parameters of the centerline extraction and 25 for full manual segmentation) and volumes of 25 CT scans were annotated by an experienced observer (ZSG). Besides, an experienced physician (DB) indicated the diameter of the PA and AA at the level of the PA bifurcation in an additional 200 CT scans. This was done in the axial view, in the slice where both the right and left pulmonary arteries appear to be of similar size at the axial view.

Parameter optimization

The centerline extraction parameters were optimized in a grid search on 10 CT scans with manually drawn centerlines. The objective was to minimize the mean centerline distance (MCD) between the manual and automatic centerlines. The radius range of the medialness filter, the cost function weights and , and the sigma for centerline smoothing were optimized. From radius range of [5 mm, 30 mm], we obtained a medialness radius range of [8 mm, 16 mm] for pulmonary arteries and [12 mm, 24 mm] for the aorta which agreed well with the radius range reported for these vessels in the literature. , Medialness filter searches for circular structures in the defined radius range. Although the method is not very sensitive to the radius range, selecting the range exactly as the average range of the vessel can result in missing vessels with abnormalities, such as aneurysm or dilatation.

From range [1, 15], the weights for pulmonary arteries and for the aorta and for both vessels were obtained, with smoothing of mm for pulmonary arteries and mm for the aorta. The extracted values for and show that medialness has a higher effect on the centerline accuracy than lumen intensity similarity since it is more accurate for localizing the vessel, especially the vessel center. Since the centerline is extracted at the center of a circular shape, the centerline can be irregular in regions where a vessel branches off, such as in the arch. A reasonable value for smooths out these jagged areas. Although we applied a grid search to fine‐tune the parameters, with our experience, the proposed method is not very sensitive to slightly changed settings of these parameters. Relatively good segmentation results can be achieved after a few trial and error runs of different parameter sets.

A five‐fold cross‐validation on 25 CT scans (independent of those used for centerline parameter tuning) was performed in which the best parameters determined from 20 CT scans were used to segment the five left out scans. In this, the best parameter set was selected as the parameter set giving the maximum average Dice similarity coefficient (DSC). The parameters to optimize were of the Gaussian Kernel used to smooth the initial segmentation, and the smoothness penalty of the optimal surface graph cuts. From the five parameter sets extracted in the cross validation, the most frequent parameter set was selected for the rest of the experiments. The obtained optimal parameter set was of 4.4 and 2.4 mm and of 32 and 40 for the pulmonary arteries and aorta, respectively.

Seed point and centerline extraction

The centerline and the seed point extraction were validated on 25 CT scans using MCD. MCD is the average symmetric Euclidean distance of all points of the automatically extracted centerline to the manual centerline.

Seed point extraction was assessed by computing the MCD between the automatic centerlines traced from manually placed seed points with those traced from automatically extracted seed points. Seed points that were placed on the border or outside the vessel resulted in a centerline with large MCD compared to the centerlines traced from manual seed points. We also used nonparametric Mann–Whitney U test to assess whether there is a significant difference between the MCD computed between the manual centerlines and the automatic centerlines traced from manually placed seed points and MCD between the manual centerlines and the automatic centerlines traced from automatically extracted seed points. Note that to prevent bias in extracting seed points of the validation set (25 CT scans), in the multi‐atlas registration method, the test scan is excluded from the atlases, and only a maximum of 24 atlases are used.

Failure in the centerline extraction in all cases (471 CT scans) could be automatically detected as follows. For the aorta, centerlines that never reached above the automatically extracted PA bifurcation point, or that went through the coarse initial PA segmentation, were considered as failed extractions.

For the pulmonary arteries, an additional centerline for only the pulmonary trunk was traced between the seed points at the PA trunk and the PA bifurcation point. Subsequently, the MCD between this centerline and the main centerlines was extracted and the centerlines that had a large MCD, or the ones which went through the coarse initial aorta segmentation, were considered as failed extractions.

Segmentation

The segmentation accuracy was assessed on 25 CT scans by comparing it with manual segmentations. The DSC was computed to assess the degree of spatial overlap of the automatic segmentation with the manual segmentation. The mean symmetric surface distance (MSD) was computed in millimeters between the manual and automatic segmentation surfaces.

For a larger scale, qualitative assessment, the pulmonary arteries and aorta were segmented on 436 additional CT scans, which had no manual annotations of the full volume. The quality of the segmentations was inspected visually by overlaying the segmentation of each vessel on the original CT scan. In the regions with boundary inaccuracy, the maximum distance between the segmentation surface and the vessel boundary was manually measured in the axial plane by ZSG. Then, based on visual inspection, the segmentations were separated into three groups: high‐quality segmentation (without significant boundary inaccuracy); segmentation with minor errors (max 3 mm SD); and segmentation requiring correction (more than 3 mm SD).

Repeatability, diameters, and PA:AA ratio

The automatic diameter of the PA and AA at the level of the PA bifurcation as well as the ratio were computed on 436 CT scans. Out of 436 CT scans, the accuracy of the in‐slice diameters and the ratio were assessed on 200 CT scans by comparing them with the in‐slice manual diameters and manual ratio.

In the remaining 236 CT scans, the average volumetric diameters and the ratio extracted from the 10 mm segment around the landmark level were computed. These 118 short‐term repeat scan pairs (236 scans) were used for assessing the repeatability of the method. Changes of the main PA diameters and the aortic diameters within the 3‐month period are expected to be negligible since the changes in the main PA diameter within 8 months are mm and the annual change in aortic diameters is 0.1–0.2 mm.

The diameters were expressed as mean standard deviation (range). The repeatability and the agreement between the manual and automatic diameter and the ratio were assessed by the intraclass correlation (ICC) [95% confidence interval], based on a single‐rating, absolute‐agreement, two‐way mixed‐effects model. The quality of the vessel segmentation on all 436 CT scans was assessed quantitatively with the visual inspection as explained in Section 3.4.

RESULTS

MCDs between the manual and automatic centerlines, for both manual and automatically extracted seed points, are shown in Table 2. The centerlines were always extracted inside the vessel and close to the vessel center. The average distance was less than 0.5 mm between the automatic centerlines traced from the automatically extracted seed points and the ones traced from the manually placed seed points. The nonsignificant difference in MCD to the manual centerlines between these two sets infers the reliability of the seed point extraction method.

TABLE 2

MCD (mm) between automatic and manual centerlines, in 25 CT scans for the left (LPA) and right pulmonary arteries (RPA) starting from pulmonary trunk (PAT), and for the aorta from both manual and automatically extracted seed points

	PAT to RPA	PAT to LPA	Aorta
Manual SP	2.14.±.0.63	2.38.±.0.66	1.59.±.0.51
Automatic SP	2.46.±.0.81	2.77.±.1.41	1.54.±.0.33
p‐value	0.06	0.27	0.35

Abbreviation: SP, seed point.

Note: p‐values are compared using nonparametric Mann–Whitney U test.

Out of all 461 CT scans used in this study, centerline extraction failed in only 17 cases (17 out of 922 vessels including the PA and aorta [1.8%]) and all failure cases were detected automatically as described in Section 3.3. In five cases, the centerline failure was due to wrong seed point extraction in which the seed points were extracted in the border or outside of the vessel. In four cases, the aortic centerline made a shortcut through the PA, and in the remaining eight cases, the PA centerline was either at the border of the vessel or made a shortcut through the background. These failures were mainly due to unclear vessel borders where the PA is adjacent to the aorta or other structures. Scans with failed centerlines were excluded from subsequent analysis.

The vessel segmentation of the 436 CT scans with no manual annotation was quantitatively assessed with visual inspection. Table 3 shows the assessed quality of the aorta and pulmonary arteries segmentation of the 419 CT scans with no failure in the centerline extraction.

TABLE 3

Qualitative assessment of the PA and aorta segmentation of 419 CT with visual inspection

Visual inspection	Aorta	PA
High quality	394 (94%)	387 (92%)
Minor error	18 (4%)	24 (6%)
Correction required	7 (2%)	8 (2%)

Note: Segmentations are with no obvious error (high quality), with max 3 mm SD (minor error), or with more than 3 mm SD (correction required).

Box plots of segmentation DSC and MSD, for the pulmonary arteries and aorta obtained in five‐fold cross‐validation on 25 CT scans with manual segmentation, are shown in Figure 3, with an average and average mm for pulmonary arteries and and mm for the aorta. Figure 5 illustrates three examples of segmentation results overlaid with manual annotations.

FIGURE 3

Box plots of the DSC and MSD between the manual and automatic PA and aorta segmentation. The plot shows the median (green), interquartile range (boxes), 99.3% coverage of the data (whiskers), and the outliers (+ symbol in red)

FIGURE 5

Examples of the obtained segmentations. First row is a good segmentation with DSC = 0.95 for PA and DSC = 0.97 for the aorta. The next two rows are the two cases with the largest segmentation errors (lowest DSC) for the PA and aorta, respectively, with DSC = 0.90, DSC = 0.96 for row 2, and DSC = 0.92, DSC = 0.95 for row 3. Column (a) shows an axial slice close to the heart and column (c) shows a slice at pulmonary artery bifurcation level. In columns (b) and (d), the manual PA (green) and aorta (yellow) are overlaid with the automatic segmentation of the PA (blue) and aorta (red). The 3D segmentation of the aorta (e) and PA (f) is shown together with the 3D view of the CT scan in (g)

From 200 CT scans used to compare manual and automatic in‐slice diameters, 11 scans failed the centerline or seed point extraction. High agreement between the diameters was obtained on the remaining 189 CT scans. The scatter plots of these diameters are illustrated in the first row of Figure 4. The first three rows in Table 4 present the diameters, diameter difference, and ICC between the automatic and manual measurements.

FIGURE 4

Scatter plots and ICC of the manual and automatic diameters measured at the axial slice at the level of PA bifurcation (first row), and of the automated average volumetric diameter estimates in scan and short‐term rescans (second row) for, from left to right, the PA, aorta, and PA:AA ratio

TABLE 4

Diameters, diameter difference, and ICC between manual and automated measurements (first three rows) and between the short‐term scan–rescan pairs (last three rows), of the pulmonary artery, ascending aorta, and PA:AA ratio

Method		Diameter (mm)		Difference	ICC [95% CI]
		Manual	Automatic
2D Axial n = 189	PA	26.39 ± 3.28 (18.32, 36.19)	25.79 ± 3.53 (18.06, 36.72)	0.60 ± 1.36	0.92 [0.89−0.94]
	Aorta	34.83 ± 3.63 (25.55, 44.53)	34.39 ± 3.74 (25.21, 43.81)	0.44 ± 0.96	0.97 [0.95−0.97]
	PA:AA	0.76 ± 0.09 (0.51, 1.01)	0.76 ± 0.10 (0.50, 1.09)	0.01 ± 0.04	0.90 [0.87−0.93]
		Scan (baseline)	Rescan
3D Volume n = 112	PA	28.24 ± 3.37 (21.58, 38.14)	28.40 ± 3.41 (21.66, 38.48)	0.16 ± 1.45	0.89 [0.87−0.93]
	Aorta	34.72 ± 3.34 (27.65, 44.22)	34.76 ± 3.35 (27.28, 45.37)	0.04 ± 0.97	0.95 [0.93−0.97]
	PA:AA	0.82 ± 0.10 (0.60, 1.12)	0.82 ± 0.09 (0.61, 1.14)	0.01 ± 0.05	0.86 [0.82−0.90]

Abbreviation: CI, confidence interval.

From 118 subjects used to assess the repeatability, six subjects had failed centerline or seed point extraction. For the remaining 112 subjects (224 CT scans), the automatic average diameter extraction from 3D volume showed a high correlation between scan and rescan of each subject with an ICC of 0.89, 0.95, and 0.86 for the PA, aorta, and the ratio, respectively. The diameters, diameter difference, and ICC between the scan–rescan pairs of 112 subjects (224 CT scans) are shown in the last three rows of Table 4 and the scatter plots are illustrated in the second row of Figure 4.

DISCUSSION

In this study, we presented a fully automatic segmentation and diameter measurement method to segment the pulmonary arteries and aorta and to measure their diameters in non‐ECG‐gated, noncontrast CT scans. Automatic extraction of the level of the PA bifurcation allowed automatic measurement of the ratio. We verified the quality of the segmentations on 25 CT scans by comparing them with full 3D manual annotations. The segmentation algorithm performed well with an average Dice overlap of 0.93 and 0.96 and mean surface distance of 0.62 and 0.43 mm, less than the in‐plane voxel size, for the pulmonary arteries and aorta, respectively.

Visual inspection indicated that 92% of pulmonary arteries and 94% of the aorta were segmented with high quality and with no obvious error. In the remaining almost 6–8%, segmentation errors mainly occurred in regions close to the heart. In this region, the PA trunk and the aortic root are adjacent and have similar intensity. This makes the vessel boundaries unclear, which results in segmentation errors. Furthermore, motion artifacts caused by the motion of the heart during the cardiac cycle increase the ambiguity of the vessel boundaries and make the segmentation difficult, even for experienced radiologists. Figure 5 shows the three cases, where the cases at rows 2 and 3 have the largest segmentation errors. Segmentation errors with respect to the manual segmentation are visible in an axial slice close to the heart (column b), whereas the segmentation has high accuracy at the level of the pulmonary bifurcation (column d).

Table 5 presents comparative segmentation results on the methods proposed in the literature to segment the PA and aorta on noncontrast CT scans. Compared to our previous work, in this study, besides segmenting the pulmonary arteries, we evaluate the method on a larger region for the aorta, from the diaphragm level at the descending aorta to the aortic root. Although results are not directly comparable due to differences in datasets, our method achieved higher Dice overlap and lower mean surface distance than was reported for all previous methods for both vessels. The advantage of the presented optimal surface graph‐cut method is incorporating the topology constraints into the graph structure and defining “smoothness‐penalty” edges to encourage a smooth segmentation similar in shape to the initialization. These constraints and penalties restrict the final segmentation to plausible vessel surfaces. In contrast, since in the deep learning‐based methods, , no prior information on the geometry of the vessels and no explicit shape constraints have been incorporated to the network, in regions with very unclear boundaries, unrealistic/incorrect surfaces can be produced.

TABLE 5

Segmentation results compared with literature

	Method	DSC	MSD
PA	Xie et al. 33	0.88	–
	Haq et al. 30	0.80 ≤ DSC ≤ 0.91	‐‐
	Presented method	0.94 ± 0.02	0.62 ± 0.33
Aorta	Išgum et al. 20	0.87 ± 0.03	–
	Avila‐Montes et al. 24	0.88 ± 0.05	–
	Dasgupta et al. 25	0.88 ± 0.06	–
	Tahoces et al. 26	0.95	–
	Xie et al. 27	0.93 ±0.01	1.39 ± 0.19
	Kurugol et al. 28	0.92 ± 0.01	0.62 ± 0.09
	Trullo et al. 29	0.89 ± 0.04	–
	Noothout et al. 22	0.91 ± 0.04	1.32 ± 0.85
	Haq et al. 30	0.75 ≤ DSC ≤ 0.94	–
	Lartaud et al. 21	0.92 ± 0.02	1.25 ± 0.78
	He et al. 23	0.95	–
	Sedghi Gamechi et al. 39	0.95 ± 0.01	0.56 ± 0.08
	Presented method	0.96 ± 0.01	0.43 ± 0.07

The running time of our presented method is about 60 min, where the seed point extraction process takes about 20 min. Therefore, although our presented method's high accuracy and reliability make it suitable for clinical study, the relatively long processing time could limit the uptake of the method in clinical practice. However, this processing is still substantially faster than full manual annotation and may reduce human errors without requiring any manual interaction.

We achieved a high agreement between the automatic and manual diameters for the PA, aorta, and ratio in axial slices, with an ICC of 0.92, 0.97, and 0.90, respectively. This is higher than the interobserver agreement reported by Terzikhan et al. of 0.91 and 0.94 for the PA and aorta diameter. Table 4 indicates that the manual diameters are slightly (on average 0.5 mm) larger than automatic diameters for both PA and aorta. This may be explained by the fact that observers do not necessarily choose the point of maximum intensity gradient as the boundary. The difference between the automated method and observer annotations is similar to the interobserver bias of 0.4 mm reported by Tonelli et al. for PA diameter on noncontrast CT scans, and interobserver bias of 0.5 mm reported by Quint et al. for mid‐AA diameter on CTA scans.

The manual diameter measurements, such as presented in Refs. 5, 8, 9, and 11, are subjective based on slice location. Measurements that are further away from the bifurcation show smaller diameters than those close to the pulmonary bifurcation. Linguraru et al. show that a small shift along the PA centerline can lead to diameter changes up to 20%. Moreover, the orientation of the vessels with respect to the patient and with respect to the axial plane may vary, leading to variability in axial diameter measurements. Also, determining the bifurcation level is difficult and is prone to variability. The 3D volumetric average diameter measurement in segments perpendicular to the vessel centerlines as proposed in section 2.5 is less subjective and is a more robust and reproducible technique than diameter measurement in 2D axial slices, which potentially decreases discordant measurements. The diameters extracted by the presented 3D method showed high scan‐rescan repeatability (Figure 4). This improved robustness and reproducibility may lead to better association with pulmonary hypertension; however, this requires further investigation.

Considering that lung cancer screening with noncontrast CT is becoming more common, our method can be used to screen for (mild) PA and aortic dilatation as biomarkers of cardiovascular disease in the same populations. A limitation of this study is that the method is evaluated on data from a single scan protocol and a relatively healthy screening population. For application in data from very different scan protocols, parameters may need to be adjusted.

CONCLUSION

A fully automatic method is presented to segment the pulmonary arteries and aorta and quantify their diameters on non‐ECG‐gated, noncontrast CT scans. Qualitative and quantitative analysis demonstrates that our method provides robust, accurate, and reproducible measurements of the PA and aorta diameters and the ratio. Automatically extracting a full 3D shape and size of the vessel, the vessel diameters, and biomarkers with high accuracy, in non‐ECG‐gated noncontrast CT scans such as are acquired in lung cancer screening, can provide important prognostic information and enable the early‐stage diagnosis of cardiovascular disease and provide factors for risk assessment in patients with COPD.

CONFLICT OF INTEREST

The authors have no conflict to disclose.