LTSP: long-term slice propagation for accurate airway segmentation.

PURPOSE: Bronchoscopic intervention is a widely used clinical technique for pulmonary diseases, which requires an accurate and topological complete airway map for its localization and guidance. The airway map could be extracted from chest computed tomography (CT) scans automatically by airway segmentation methods. Due to the complex tree-like structure of the airway, preserving its topology completeness while maintaining the segmentation accuracy is a challenging task.
METHODS: In this paper, a long-term slice propagation (LTSP) method is proposed for accurate airway segmentation from pathological CT scans. We also design a two-stage end-to-end segmentation framework utilizing the LTSP method in the decoding process. Stage 1 is used to generate a coarse feature map by an encoder-decoder architecture. Stage 2 is to adopt the proposed LTSP method for exploiting the continuity information and enhancing the weak airway features in the coarse feature map. The final segmentation result is predicted from the refined feature map.
RESULTS: Extensive experiments were conducted to evaluate the performance of the proposed method on 70 clinical CT scans. The results demonstrate the considerable improvements of the proposed method compared to some state-of-the-art methods as most breakages are eliminated and more tiny bronchi are detected. The ablation studies further confirm the effectiveness of the constituents of the proposed method and the efficacy of the framework design.
CONCLUSION: Slice continuity information is beneficial to accurate airway segmentation. Furthermore, by propagating the long-term slice feature, the airway topology connectivity is preserved with overall segmentation accuracy maintained.

Introduction

Bronchoscopic intervention [1-3] is widely used in clinical practice since it could achieve minimally invasive access techniques for pulmonary diseases. To improve the efficacy of the bronchoscopic intervention, an accurate 3D airway map extracted from chest computed tomography (CT) is required, which is essential to the endoscopic tips localization and intraoperative guidance. However, due to the complex tree-like structure, manual segmentation of the airway from chest CT scans is time-consuming and requires expert knowledge. Therefore, to relieve the burden of clinicians, automatic airway segmentation methods are proposed to extract airways accurately.

Traditional airway segmentation methods are always designed based on region growing [4, 5] and morphology operation [6], which are very sensitive to the manually designed features. Thus, their performance will be degraded once noises appear. Besides, owing to the similarity between bronchi and surrounding human tissues, these methods will cause severe leakages and breakages. As shown in Fig. 1, in bronchoscopic intervention practice, leakages in the airway segmentation will lead to misleading destinations while breakages will cause interrupted trajectories. Recently, convolutional neural networks (CNNs)-based methods [7-11] are increasingly developed for airway extraction. U-Net [12] or 3D U-Net [13] is widely used as the backbone to acquire a coarse airway segmentation result, which is then refined by utilizing prior knowledge like neighborhood connectivity [14] or global information [15].

Fig. 1

The impact of leakages and breakages in the airway for bronchoscopic intervention. a describes the principle of bronchoscopic intervention referred to [3]. b demonstrates that airway leakages will provide an incorrect destination and mislead the trajectory. c illustrates that airway breakages will interrupt the planning trajectory

However, due to the complex tree-like structure of the airway, acquiring accurate and fine-grained segmentation results is difficult. There are several challenges to be solved in accurate airway segmentation tasks. First, the intensity distribution is different in the main trachea region and bronchi region. In the main trachea region, the intensity contrast between the airway lumen and the wall is distinct and the features are easy to learn, while the ambiguity of intensity contrast in the bronchi region is hard to identify for the CNN models. Second, the encoder–decoder architecture is widely used in the CNN models like U-Net, where several pooling operations are utilized. Since some bronchi only have a diameter of 2-5 voxels, these features will be vanished by pooling operations and are hard to reconstruct in the decoder. Furthermore, these two challenges will result in weak airway features in the decoding stage and affect the accuracy of airway segmentation. Since adjacent slices in chest CT scans have similar airway shape and lumen position, one intuition for weak feature enhancement is to use this airway continuity prior knowledge. The prior knowledge could be achieved by applying slice feature propagation, so as to strengthen the weak feature and impose the continuity constraints.

To address the above problems and preserve the topological connectivity of the airway, we propose a two-stage end-to-end framework for accurate airway segmentation using long-term slice propagation (LTSP). In the proposed method, the slice relationship is considered in the decoding stage of the framework. The slice feature is transferred in the proper direction to fully recover the continuity information destroyed in the encoder. Although slice propagation could alleviate the breakage phenomenon by enhancing the airway features, there remains a problem that slice features are hard to transfer among too long distances. In our method, to solve the above gradient vanish problem, we further design the LTSP cell in the decoding stage where continuity information is extracted and densely propagates to other adjacent slices. Furthermore, we compare the proposed airway segmentation method to some state-of-the-art methods in 70 clinical CT scans. Extensive experiments show that our method achieves superior performance in extracting topological complete airway while maintaining the competitive overall segmentation accuracy.

Methodology

The overview of the proposed airway segmentation method is illustrated in Fig. 2, which is a two-stage end-to-end framework. The first stage is used to acquire a coarse feature map from the input cube by utilizing an encoder–decoder architecture. The second stage is to predict a refined and well-connected airway segmentation result by: (1) propagating slice information of the coarse feature map using the long-term slice propagation (LTSP) method, (2) decoding the final airway segmentation result from the refined feature map with skip connection. The whole pipeline processes the 3D input cube as a whole volume except the LTSP module processes it sequentially.

Fig. 2

Illustration of the proposed airway segmentation framework. The channel number is denoted below each feature map. In the first stage, a coarse feature map is extracted from 3D CT cropped cube by an encoder and a coarse decoder. In the second stage, each slice feature in the coarse feature map is propagated by passing through the LTSP cell. The results are stacked together and then used to predict the refined segmentation result by a fine decoder

Stage 1: coarse feature map generation using encoder–decoder architecture

Stage 1 in the framework is to provide a coarse feature map by using an encoder and a coarse decoder. The feature map is extracted for subsequent slice propagation and airway refinement.

In this stage, we employ an encoder–decoder architecture to extract the airway feature and then recover from the deep information. A part of 3D U-Net is used as the backbone, which contains three down-samplings and two up-samplings with skip connections. In each down-sampling module, two convolution layers (Conv) with batch normalization (BN) and rectified linear unit (ReLU) are followed by a max-pooling layer. In each up-sampling module, the up-sampled feature map and corresponding feature map in the encoder are concatenated together and then passed to a convolution layer with BN and ReLU.

Given a 3D CT cube input with size of , a coarse feature map with size of is generated by stage 1 aswhere denotes the feature extraction process in stage 1.

Stage 2: segmentation refinement using long-term slice propagation

Stage 2 in the framework is used to densely propagate the slice features sequentially in the coarse feature map extracted in stage 1. The proposed LTSP method is designed to transfer the long-term slice features for the refined feature map generation. The final airway segmentation result is then predicted from the refined feature map by a fine decoder.

Long-term slice propagation

Taking the airway’s continuity prior knowledge into consideration, we utilize the slice propagation method into the decoder, which aims to enhance the weak airway features in each slice by integrating the adjacent slice features. Additionally, we improve the previous spatial CNN [16] method and propose a long-term slice propagation (LTSP) method for effective feature transferring. For the place of LTSP module, as skip connection enables features in the decoder to integrate with features in the encoder, placing the LTSP module in the decoder could more effectively impose continuity constraints. Therefore, the airway features are more efficiently recovered by the guidance of features in the encoder.

The original method proposed in [16] uses the spatial CNN to propagate the neighboring message in 2D images. However, as it only transfers the feature from a slice to its nearest neighbor, the propagating airway feature will be rapidly dismissed and the continuity information will be lost. To improve the propagation efficiency and expand the receptive field of slice propagation, we propose the long-term slice propagation (LTSP) method.

Figure 3 gives the comparison of spatial CNN method and LTSP method. As illustrated in Fig. 3a, in spatial CNN, one slice only receives the feature from its upper slice by directly utilizing convolution operation. For the complicated tree-like airway structure, there are several limitations in spatial CNN: (1) insufficient feature propagation. Some breakages in the airway are dependent on the feature far from it, which needs to use features in several slices to recover it. As one slice only can receive the feature from its nearest neighbor, the information is not sufficient for its recovery. (2) Short propagation distance. Spatial CNN only uses convolution to extract the information from a neighbor slice, leading to the feature disappearance when the propagation distance is too long. Besides, a short propagation distance will result in a narrow receptive field and great propagation deficiency.

Fig. 3

Comparison of spatial CNN method and long-term slice propagation (LTSP) method. a gives the principle of spatial CNN where convolution operation is used to propagate the slice features in coarse feature map slice-by-slice. b illustrates the principle of long-term slice propagation. Each slice in the coarse feature map is passed through the LTSP cell to acquire a continuity information map which is then fused with the original feature by slice propagation. c demonstrates the mechanism of the designed LTSP cell. It receives the current slice feature , previous output and previous cell state to acquire the current output and cell state by using the updating rules

To address the problems, we propose the LTSP method to develop the effectiveness of slice propagation. As illustrated in Fig. 3b, the LTSP method consists of two parts: information extraction and feature propagation.

In the information extraction part, the LTSP cell inspired by the long short-term memory (LSTM) [17] is designed to extract airway continuity information. Figure 3c gives the propagating mechanism of the designed LTSP cell. Given a coarse feature map with the size of , each split slice with the size of is acquired byTogether with the previous LTSP cell output and previous cell state (with each size of ), the forward calculation includes the following three steps. First, the stacked gates G with the size of is acquired by using 2D convolution operation to double the channel number of the concatenation of and , where G is computed byThen, the convolution result G is chunked into , , and , with each size of . Pass them into different activation function to calculate forget gate , new cell state , input gate , and output gate , respectively. The calculation process could be described asFinally, previous cell state and new cell state are used to calculate the updated output and final cell state byIn the feature propagation part, the original slice is updated to by propagating information maps and extracted by LTSP cells. The forward slice propagation could be described aswhere and denote the convolution weights for 1-distance feature and 2-distance feature, respectively. And denotes the nonlinear activation function like ReLU.

In the end, all updated slices are stacked into a new refined feature map with the size of by

Segmentation prediction and optimization

To acquire the airway segmentation from refined feature map , a fine decoder with skip connection is used to get the probability map of the airway . An argmax operation is then utilized to predict the final airway segmentation result . The prediction process could be described aswhere denotes the up-sample operation in fine decoder.

In the training process, soft dice loss [18] is used for the airway segmentation tasks. Given prediction p(x) and corresponding binary label y(x) for each voxel x in segmentation result volume X, the segmentation loss could be calculated bywhere smoothing parameter is used to avoid division by zero.

Experiments and results

We evaluate the proposed method in LIDC dataset [14], where 50 chest CT scans are randomly chosen for training and the remaining 20 scans are used for testing. Furthermore, ablation studies are conducted to confirm the effectiveness of our method.

Datasets and implementation details

The experiment dataset contains 70 clinical chest CT scans, where the pixel spatial resolution ranges from 0.5 to 0.781mm and slice thickness varies from 0.45 to 1.0mm. The model is trained on 50 randomly chosen scans and tested on the remaining 20 scans.

To improve the model’s generalization ability, the HU value of each scan is truncated into [−1000, 600] and linearly mapped into [0, 1]. In the training phase, each CT cube of size is cropped from original CT scans as the input of the model. Besides, data augmentation is performed to each cropped cube including horizontal flipping and slight rotation. For training strategy, we adopt Adam optimizer (, ) with learning rate set as 0.002. Our method is implemented in PyTorch 1.7 with NVIDIA GeForce RTX 3090. The training process converges in 30 epochs.

Additionally, we utilize the center crop method instead of the random crop method to acquire more representative cropped cubes as the model’s inputs. The center crop method is used to guarantee that the cropped cubes are able to contain sufficient airway features. Given a CT scan with a size of , we first search for the minimum and maximum index for airway region in x, y, and z directions. We denote these index pairs as and randomly choose a voxel from these spans as the cropped center. The input cube and corresponding binary label with a size of are cropped from original CT scans and ground truth maps based on the cropped center. Then, the cropped cubes are fed into the model for training. Furthermore, we utilize the sliding window technique in the testing process to acquire each cube’s prediction and combine them to form the whole airway segmentation result.

Evaluation metrics and results

To evaluate our method, we adopt four metrics to assess the topological completeness and segmentation accuracy of airway prediction results: (1) branches detected (BD), (2) tree-length detected (TD), (3) Dice coefficient (DSC), and (4) false positive rate (FPR). The definition of four metrics could be referred to [19].

In our method, 3D U-Net is adopted to be the backbone of the segmentation framework. Therefore, we compare our method with original 3D U-Net [13] and other related state-of-the-art methods like Wang et al. [20] and Juarez et al. [15]. Wang et al. [20] develop a radial distance loss for detecting more tiny airway tubular structures. Juarez et al. [15] introduce graph neural network (GNN) module into the deepest level of 3D U-Net to improve airway connectivity. These methods are implemented by ourselves and fine-tuned on the experiment dataset. Table 1 gives the comparison results which show that our method achieves the best performance in BD and TD with compelling DSC and FPR. Compared to others, our method increases the BD and TD by over 3 % and 2 %, respectively. Since these two metrics could reflect the topological completeness of airway segmentation, the comparison results also demonstrate that our method outperforms the others in detecting more small bronchi and improving the connectivity of airway segmentation.

Table 1

Results (%) of the proposed framework compared to state-of-the-art methods (Mean±Standard deviation)

Method	BD	TD	DSC	FPR
3D U-Net [13]	86.06±11.72	83.10±9.62	93.35±1.74	0.02±0.01
Wang et al. [20]	87.97±9.49	84.61±8.54	93.12±2.00	0.02±0.01
Juarez et al. [15]	87.40±10.15	84.60±8.94	93.51±1.57	0.02±0.02
Our proposed	90.83±9.38	87.59±8.71	92.95±1.61	0.03±0.01

The figure in bold style indicates the best performance in each metric

Qualitative comparison of airway segmentation shown in Fig. 4 also demonstrates the effectiveness of our method. Compared to other methods, more small branches and tiny bronchi are reconstructed in our method, which results in the great improvement of the airway topological connectivity. The preservation of airway topological completeness is significant for clinical practice like bronchoscopic intervention and surgical navigation. In Fig. 4b, a great number of breakages (false negatives voxels in green color) are eliminated in our method while maintaining a high overall segmentation accuracy.

Fig. 4

Rendering of airway segmentation results. a and b give the comparison of different methods in an easy case and a hard case, respectively. The true positive voxels are shown in red color, while the false negative voxels are shown in green color

Moreover, to assess the generalizability of the proposed framework on hard cases with abnormality, we train our model in the LIDC dataset and directly test it on the coronavirus 2019 (COVID-19) dataset [10], which contains 58 CT scans with pulmonary diseases. Figure 5 gives the comparison of airway segmentation results with 3D U-Net baseline and the proposed framework on a specific severe case with a huge abnormality. Qualitative results in 3D volume illustrate that the proposed method could detect more tiny bronchi and enhance the overall topological connectivity. Comparisons on 2D planes further validate the effectiveness of the proposed method in detecting airway lumen influenced by abnormality.

Fig. 5

The qualitative comparison in 3D volume and 2D planes of baseline and the proposed LTSP method on a severe case in the COVID-19 dataset [10]. The true positive voxels are shown in red color, while the false negative voxels are shown in green color

Ablation study

We also conduct ablation studies to further investigate the effect of each component in our method and the impact of placing the LTSP module in the different stages of the framework. Since we use dense structures in the LTSP method, the density of this structure should be considered in our framework. Therefore, experiment comparisons are made by utilizing modules with different propagation distances in a single pass. To validate the effectiveness of the LTSP cell in slice propagation, experiments with or without LTSP cells are also conducted. Furthermore, we change the place of the designed LTSP module from decoding stage into encoding stage and bottleneck to compare the difference.

As shown in Table 2, when placed in the decoding stage of the framework, the LTSP cells and more densely connected structures both can achieve higher BD and TD with DSC maintained. For the comparison of propagation distance, the two-slice module outperforms the no propagation module and one-slice module since it enables the slice feature to propagate to more slices directly. Furthermore, experiment results show that the designed LTSP cell could effectively improve the segmentation accuracy and topological completeness by extracting more valuable continuity information.

Table 2

Comparisons (%) of propagation density and LTSP cells on the testing set when the LTSP module is placed on the decoding stage of the proposed framework (mean±standard deviation)

Module configuration	BD	TD	DSC	FPR
No propagation module	86.06±11.72	83.10±9.62	93.35±1.74	0.02±0.01
One-slice module	87.97±9.49	84.61±8.54	93.12±2.00	0.02±0.01
One-slice + LTSP cells	90.11±10.28	87.20±8.96	93.15±1.79	0.02±0.02
Two-slices module	89.22±9.41	86.39±7.84	92.42±2.79	0.03±0.02
Two-slices + LTSP cells	90.83±9.38	87.59±8.71	92.95±1.61	0.03±0.02

The figure in bold style indicates the best performance in each metric

Qualitative results in Fig. 6 also demonstrate the effectiveness of increasing direct propagation distance and the adding value of LTSP cells. Compared with the one-slice module, more tiny branches are recovered and some breakages are eliminated by applying the two-slice module. Besides, introducing LTSP cells into the one-slice module is also beneficial to detecting more essential branches and slight bronchi. The development of topological completeness is significant in clinical applications like bronchoscopic intervention.

Fig. 6

Comparisons of airway segmentation results on a specific case under different LTSP module configurations. The true positive voxels are shown in red color, while the false negative voxels are shown in green color

Moreover, comparisons are also made by placing the LTSP module in the encoding stage, bottleneck, and decoding stage in the framework with the same configuration. Table 3 gives the results which demonstrate the effectiveness of placing the LTSP module in the decoding stage of the framework. The quantitative results also support the theoretical knowledge that the LTSP method could use more robust information to enhance the slice feature since the feature map in the decoding stage is integrated with the feature in the encoding stage by skip connection.

Table 3

Comparisons (%) of placing LTSP module in the different stages in the framework on the testing set (mean±standard deviation)

Stage	BD	TD	DSC	FPR
Encoding stage	87.00±12.40	83.92±10.44	92.57±2.67	0.02±0.02
Bottleneck	87.48±12.52	84.55±10.18	93.34±1.90	0.02±0.01
Decoding stage	90.83±9.38	87.59±8.71	92.95±1.61	0.03±0.02

The figure in bold style indicates the best performance in each metric

Conclusion

This paper proposed an effective long-term slice propagation (LTSP) method for accurate airway segmentation. Focusing on alleviating breakage phenomenon and improving topological connectivity, LTSP cell and dense slice propagation were designed to fully exploit the slice continuity relationship. Extensive experiments showed that the proposed method outperformed some state-of-the-art methods by detecting more tiny bronchi and reconstructing essential branches, which further validates the effectiveness of its constituents and framework design. The proposed method was also assessed in severe cases with abnormality to validate its generalizability. The proposed method is beneficial to extracting topological complete airway segmentation for bronchoscopic intervention.