Automatic three-dimensional segmentation of endoscopic airway OCT images

: Automatic delineation and segmentation of airway structures from endoscopic optical coherence tomography (OCT) images improve image analysis efficiency and thus has been of particular interest. Conventional two-dimensional automatic segmentation methods, such as the dynamic programming approach, ensures the edge-continuity in the xz-direction (intra-B-scan), but fails to preserve the surface-continuity when concerning the y-direction (inter-B-scan). To solve this, we present a novel automatic three-dimensional (3D) airway segmentation strategy. Our segmentation scheme includes an artifact-oriented pre-processing pipeline and a modified 3D optimal graph search algorithm incorporating adaptive tissue-curvature adjustment. The proposed algorithm is tested on endoscopic airway OCT image data sets acquired by different swept-source OCT platforms, and on different animal and human models. With our method, the results show continuous surface segmentation performance, which is both robust and accurate.

submucosa-cartilage interface through either using a drawing tablet, or using a mouse to label landmark points followed by curve fitting. Since a complete airway OCT data set may have thousands of images, this manual labeling procedure is usually time-consuming. It also suffers from inter-observer variation because manual operations are subjective.
Various approaches have been proposed to tackle this problem by using automatic computer algorithms, such as the clustering algorithm [23], morphological operations [24], and the shortest path methods [21] and so on. Recently, we proposed an automatic segmentation method based on the dynamic programming (DP) algorithm to delineate the airway wall and to measure the thickness of the mucosa and submucosa layers [25]. The DP algorithm was also successfully applied to endoscopic esophageal OCT, as presented in [26,27]. However, all these methods are two-dimensional (2D) approaches. When segmenting a 3D image data set, these approaches process each image slice separately, yielding a set of individual contours rather than a complete tissue surface. The 2D methods often result in discontinuous artifacts in the y-direction (inter-B-scan direction) since the algorithm does not consider inter-slice continuity. Moreover, due to the same reason, 2D methods will be easily affected by other artifacts in the images, such as overly-saturated pixels, speckle noise, and mirror images, etc.. In these cases, manual intervention and pre-conditioning of the image data is often required.
Three-dimensional segmentation methods that consider the inter-slice connection may be used to tackle these problems. In fact, 3D surface detection algorithms have already been proposed for the segmentation of retinal OCT images, such as the 3D-expansion DP method [28], and the 3D optimal graph search algorithm [29][30][31]. However, unlike retinal OCT, which is based on raster scanning of the laser beam, the images in endoscopic OCT are acquired by rotating and pulling back of a miniature fiber-optic probe, and thus have some distinct features. The first one is various imaging artifacts ( Fig. 1(a) and 1(b)), including the plastic sheath used to protect the imaging probe and mirror images induced from active phase modulation, etc.. Second, as discussed above, the curvature of the tissue surface is usually much larger than retinal OCT or other endoscopic OCT, and it changes dramatically across the image ( Fig. 1(b)). Third, when the plastic sheath contacts the tissue surface, it will become indistinguishable ( Fig. 1(c)), further adding difficulties to separate the airway tissue. Last but not least, in order to achieve a long imaging range while preserving the pixel resolution, airway OCT images have to have more A-line pixels than conventional OCT images, leading to much higher computational demand. Apart from these features, airway OCT images also shares other universal OCT problems including speckle noise and overlysaturated A-lines due to back reflections. All these features together make it difficult to directly apply the well-developed 3D segmentation algorithms. Fig. 1. Typical endoscopic airway OCT images. Various artifacts were presented: (a) normal case with flat tissue surface (the imaging probe is in the center of the airway), noise and saturated A-lines are presented; (b) curvature of the airway wall changes dramatically when the probe is not in the center of the airway; (c) edge blurring is induced when the probe sheath is in touch of the airway tissue (inset shows the anatomy of the airway wall).
Addressing these challenges, we present a 3D processing scheme tailored for the segmentation of endoscopic airway OCT images. In the core of our method is the 3D optimal graph search (GS) algorithm originally presented in [32]. To accommodate the features in the airway OCT images, the algorithm is modified by incorporating the prior information derived from the slope of the tissue surface, and thus achieves a more robust performance. The proposed scheme also features an artifact-oriented pre-processing step, which removes the imaging artifacts such as the speckle noise and the plastic sheath, and at the same time obtains the surface slope of the airway lumen. This slope information is further used as a constraint for the subsequent graph search problem. The airway lumen, the mucosa/submucosa boundary, and the submucosa/cartilage boundary are detected by applying the slopeconstrained 3D graph search algorithm. The proposed method was tested on airway OCT data sets of animals and human obtained from two different OCT systems. Qualitative and quantitative comparisons of the new method to the 2D DP approach, the conventional 3D GS approach and the manual-labeled method were carried out. The results show that our algorithm has a robust and accurate performance, with the continuity of the tissue surfaces well preserved.

Algorithm overview
A flow chart outlining the main steps of the algorithm is shown in Fig. 2. The proposed segmentation method is composed of two major steps: 1) image preprocessing and 2) surface segmentation. In the pre-processing step, image noise and other artifacts are removed by performing denoising, sheath identification and removal, and binarization and region filtering. After detecting the upper boundary, the slope of the airway is obtained, and will be further used to constrain the subsequent 3D graph search problem. In the segmentation step, the airway lumen is first detected, and then used as a reference to flatten the structures under the airway surface. This flattening procedure will dramatically reduce the size of the graph search problems for segmenting the mucosa and submucosa layers, and thus will be able to largely cut down the total processing time. Finally, the segmentation results can then be 3D reconstructed for visualization and analysis.

Image pre-processing
The purposes of this step are: 1) to suppress and remove noise and other artifacts on the airway OCT images and 2) to extract the slope information of the airway structure to improve subsequent surface segmentation. The pre-processing consists of the following steps: 1) Denoising: the bilateral filter is applied to perform speckle noise suppression.
Compared to the medium filter, the bilateral filter achieves a better performance in terms of noise suppression and edge preservation.
2) Removing the sheath: the prominent plastic sheath causes strong disturbance in the detection of airway structure. In order to remove the sheath, the upper boundary and lower boundary of the plastic sheath in the OCT images are identified by applying conventional 3D graph search algorithm while limiting the search region. The sheath is then removed by cropping the image from the lower boundary of sheath.
3) Binarization and region filtering: a custom threshold value is used to transform the image to binaried image, and then the regions containing the entire airway structure is detected from the binarized image by filtering the image using a minimum area value and a maximum number of connected regions. This will filter out artifacts such as the mirror image artifacts.

4)
Detecting the upper boundary: the detected binarized airway region was considered to contain reliable airway structures, thereby the upper boundary of this region is regarded as a rough estimate of the tissue surface. The single surface 3D graph search method is carried out to detect this boundary, and the slope at each point along the boundary is calculated and saved for subsequent surface segmentation.
Notice that the detected boundary in this step is only utilized to obtain the surface slope.
The intermediate results during the above pre-processing steps are illustrated in Fig. 3. An effective pre-processing offers the slope information of the airway and deals with artifacts in the images, providing a solid basis for the accurate segmentation of the airway structure.

Slope-constrained graph search
The 3D optimal graph search algorithm, first proposed by Li, et. al. [32], is an effective 3D segmentation method designed for volumetric image data. In the original algorithm, a volumetric image of size X × Y × Z is considered as a 3D matrix I(x, y, z). A graph G = (V, E) consists of a set of vertices V and edges E is constructed according to this matrix. As illustrated in Fig. 4(a) and 4(b), each node in V corresponds to each voxel in I(x, y, z) and the edge set E consists of intra-column arcs and inter-column arcs that connect the neighboring voxels. The surface in a graph can be defined by a function z = f (x, y), where x ∈ {1, 2, .., X}, y ∈ {1, 2, .., Y}, and z ∈ {1, 2, .., Z}. Associated with each (x, y) pair is a column of voxels, and there is only one voxel, the voxel (x, y, f (x, y)), intersects with the surface. Cost values representing the similarity of the voxels belonging to a surface are assigned to the edges, and then the process of finding the optimal surface is transformed to solving a minimum s-t cut problem, which can be efficiently solved by the max-flow/min-cut algorithm [33].
During the graph construction process, two smoothness constraints ∆x and ∆y are used to limit the height variation of two neighboring surface points along the x and y direction ( Fig.  4(a)). In conventional optimal 3D GS algorithm, the smoothness constraints in the x and y directions are set to constants. However, such a setting cannot satisfy the situation where the slope of the structure is varying, such as in the airway OCT images. As shown in Fig. 4(d), setting a small smoothness constraint value leads to incorrect surface detection results in airway regions that are too steep. Conversely, by setting a large constraint value, the segmentation result will be more susceptible to noise and artifacts, as shown in Fig. 4(e). Therefore, an adaptive smoothness constraint that varies with the slope of the airway is desirable for accurate segmentation. To enforce this, in this proposed method, the slope of the upper boundary obtained from the pre-processing step is used as locally varying smoothness constraints to construct the 3D graph. In this manner, steep regions correspond to large ∆x and ∆y, while flat regions correspond to small ∆x and ∆y, such that the resulted optimal surface could better conform to the actual tissue curvature. The slope of the airway wall can be described by the tangent value at each point on the curve identified from the binarized region in the pre-processing step. Using x-direction as an example, the slope in this direction is: where x and z are the horizontal and vertical pixel coordinate values of points on the initial curve, respectively. X and Y are the total number of A-lines and B-scans respectively. Then, the smoothness constraint along the x-direction is written as, c is a constant to avoid S i being zero, and [·] is the rounding operation. S i allows for 99% of the expected changes along the z direction of two adjacent columns (assuming a Gaussian distribution). Also notice that the mucosa and the submucosa layers share similar topology (degree of bending), so the same smoothness constraint can be used for the segmentation of all these layers. Accordingly, the slope constraint in y direction, that is, ∆y between adjacent B-scans, can also be calculated following the above paradigm. Another critical step in graph construction is to design the cost c(x, y, z) according to the variations of grayscales of different layers. The cost assigned to each node should be inversely proportional to the probability for that voxel belonging to a surface, such that the optimal surface consists of the set of voxels with minimum total cost. In our method, we use the Sobel filter to obtain the cost at each node. Finally, a weight value are assigned to each node in the graph, and it is defined as: In conventional 3D graph search, the time complexity and memory occupation are proportional to the number of voxels and arcs of the graph, and thus to the image size and total number of surfaces. For endoscopic airway OCT, the image data sets often consist of a lot of B-scans, and these B-scans usually have a large number of A-line pixels (usually >1000) to preserve resolution. Therefore, the computational demand for processing the airway OCT images is high. To solve this, when detecting multiple surfaces, we simply flatten the image according to the segmented airway lumen boundary. Then the region below the lumen surface was cropped out to limit the range for segmenting the mucosa and submucosa surfaces (Fig. 5). After flattening, the mucosa/submucosa surface and the submucosa/cartilage surface are segmented separately with the slope-constraint graph search method. This simple flattening method doesn't require down-sampling of the original data set, therefore the resolution of the OCT images, and thus the segmentation accuracy, can be preserved.

Endoscopic airway OCT systems and OCT image data sets
The data sets in our experiment were acquired by two endoscopic airway OCT imaging systems, the detailed of which could be found in Ref [15]. and [17]. Figure 6 shows the schematic diagrams of these two systems. Briefly, the OCT system shown in Fig. 6(a) uses a swept-source (SS) laser with 120 nm bandwidth centered at 1310 nm, and a phase modulator that shifts the interference signal for extended imaging range [15]. The depth resolution of this system is 10 μm in tissue, and the maximum imaging range is 25 mm. For the system shown in Fig. 6(b), a swept-source VCSEL laser is utilized as the light source, achieving a very long coherence length [17]. The measured 3 dB imaging range is 30 mm and the axial resolution is 10 μm in tissue. In this paper, a sheep airway data set and a pig airway data set acquired by the first system, and an adult human airway acquired by the second system were used to test the performance of our method. These data sets had previously obtained in preclinical and clinical in vivo studies involving multiple animals and subjects. The characteristics and details of the subjects and animals in these studies have been described in [15] and [17]. Table 1 summarizes the image parameters of the testing data sets.

Segmentation of the airway lumen (single surface)
In a first step, the airway lumen was detected using the proposed method. For comparison, we implemented the dynamic programming method detailed in Ref [25]. Briefly, the DP method works by recursively finding the shortest-path (with minimal cost) across a node-weighted graph derived from the image. Because the final result of the DP method was a path, rather than a surface, the algorithm must be conducted on each slice separately. Our implementation of the DP method was written in MATLAB, and was the same as described in Ref [25]. The conventional 3D GS method was also carried out for comparison (with constant ∆x and ∆y, empirically determined). Figure 7 shows some representative segmentation results of the adult human and the sheep data sets. The insets in the figures are enlarged regions indicated by the white dashed lines. The DP method failed to detect the accurate lumen edge in various scenarios, including absence of partial structure (Fig. 7(b)), strong A-line reflection artifacts (Fig. 7(c)), and large surface curvature changes (Fig. 7(d)). Conventional 3D GS method also failed in steep tissue regions (Fig. 7(a) and 7(b), white arrows). By contrast, the proposed method performed much better in all these cases. Moreover, the proposed method could correctly identify the portion of the lumen edge when it was in contact with the sheath, as indicated by the blue arrows. This was due to the consideration of inter-slice information, and thus was more robust to local artifacts and noise. To further access the 3D segmentation results, the auto-segmentation results of the airway lumen were transformed into the circumferential images and were then imported to a commercial 3D visualization software (Amira, Thermo Fisher Scientific, MA, US) to perform 3D reconstruction. The results are shown in Fig. 8. As can be seen, compare to the DP method, the proposed method can better preserved the continuity of the surface along the pull-back direction. Incorrect detection results (white arrow) were found on the DP result because it is more vulnerable to local artifacts in the images.

Segmentation of the airway structures (multiple surfaces)
Multiple airway layers were further segmented for the sheep and the pig data sets. Figure 9 and Fig. 10 illustrate the single layer segmentation result and two layers segmentation result respectively. In Fig. 9, two consecutive images and their segmentation results are shown. As can be seen, for the DP method, large variations of the detected edge were found in the consecutive images (white arrows), while for the proposed method, the airway wall structure was accurately segmented. This was contributed by the robust 3D segmentation scheme that utilized the inter-slice information. Moreover, compared to conventional 3D GS method, the boundaries detected by the proposed method appeared closer to the airway surface in the steep regions, resembling the correct position of the airway-air interface (blue arrows). In Fig.  10, when segmenting the two layers of the airway, the proposed method also outperformed the DP and the 3D GS methods, with the mucosa/submucosa boundary (red line), and the steep edges (white arrows) identified more accurately. Fig. 9. The segmentation of the airway structure: (a) are two consecutive test images in the pig airway data set, (b) are the segmentation results for one layer (the airway structure) by using different methods. 3D reconstruction of the entire airway structure was carried out with the Amira software, and the results are shown in Fig. 11. It can be seen that the segmentation result of the 2D DP method suffered from many strange protuberances and dents across the tissue surfaces, especially on the mucosa/submucosa, and the submucosa/cartilage boundary surfaces (white arrows). On the other hand, the proposed 3D segmentation method successfully dealt with artifacts and noise in the original images, and the 3D reconstruction results appeared much smoother, more continuous and resembled the correct anatomical structure of the airway. Fig. 11. Multiple surface 3D reconstruction: (a) the original image from the pig data set; (b) and (c) are the 3D reconstruction of the entire airway segmented by the proposed method and the DP method, respectively; (d) to (f) are the 3D visualization of the airway lumen surface, the mucosa/submucosa boundary, and the submucosa/cartilage boundary detected by the proposed method, respectively; (g) to (i) are the corresponding surfaces detected by the DP method.

Quantitative comparisons
In order to quantitatively evaluate the performance of the proposed segmentation method, numerical comparisons to manual delineation, the DP algorithm, and conventional 3D GS method were carried out. In this work, manual delineation was carried out in three steps: 1) a series of landmark points on the edge (no less than 50 points per edge) were labelled by using a mouse; 2) the points were spline-interpolated to obtain an initial smooth and continuous edge; and 3) the edge was further refined by either adding landmark points or removing false points to modify the spline-interpolated curve until the observer was satisfied. Different edges were labeled separately using the same procedure. In our experiments, 20 images were randomly selected from each data sets (60 images in total) and manual segmentation of the edges was performed by two experienced airway OCT image analyst. Their labeled edges were averaged to obtain the final result.
In order to quantitatively measure the accuracy of the proposed method, three validation metrics were used: the root mean square error (RMSE), the mean absolute deviation (MAD) and the Dice's similarity coefficient (DSC). They are given by: where M is the result of manual segmentation, G is the result of a comparing approach (the proposed method, or the DP method, or the 3D GS method), n is the total number of pixels, and S a and S m are the areas of auto-and manual-segmented layer regions. Table 2 lists the obtained RMSD and MAD of the auto-segmentation methods in all three data sets. As can be seen, the proposed method outperformed the DP and conventional 3D GS methods in the segmentation of all types of edges across all data sets. The accuracy of the segmentation of the mucosa/submucosa edge was higher than that of the submucosa/cartilage boundary. This is caused by degradation of signal as light penetrates deeper. The DSC results of the sheep and pig data sets are listed in Table 3. The mean DSC of the proposed method was higher than that of the other methods, indicating more reliable layer segmentation results. Overall, comparing to the 2D DP and conventional 3D GS methods, the proposed 3D method had a more accurate and robust performance. However, the auto-segmentation accuracy on OCT images has a theoretical bound subject to the imaging physics; a further analysis for airway OCT images is desirable, as the one for retinal OCT presented in Ref [34].

Computational performance
The proposed algorithm was implemented in C + + , and was tested on a desktop computer with an Intel Core i7 processor and 24 GB of RAM. The average running time of the algorithm on the sheep and pig airway data sets were 377 s and 342 s for a single edge, 485 s and 448 s for two edges, and 630 s and 523 s for three edges, respectively. Because the DP method was implemented in MATLAB, it might not be fair to compare the computational performance of the two methods because codes executed on MATLAB are usually much slower than on C + + .

Discussion and conclusion
Automatic segmentation of the airway structures and layers could help boosting the image analysis efficiency in endoscopic airway OCT imaging. However, unlike conventional OCT systems, the processing of endoscopic OCT images is challenging mainly due to the proberotation-based scanning principle. In this work, we presented an automatic three-dimensional airway segmentation scheme based on the 3D optimal graph search algorithm. The major contributions of our proposed processing scheme are: 1) comparing to 2D segmentation methods such as the dynamic programming algorithm, our approach could preserve the 3D continuity of the airway tissue surfaces; 2) we proposed an modified version of the 3D optimal graph search algorithm that utilized the slope of the airway as a constraint to achieve better surface detection results; 3) by performing the artifact-oriented pre-processing, the algorithm also deals with various artifacts in the images, including reliable sheath removal and noise cancelling; and 4) special measures that dramatically cuts down the scale of the graph search problem are taken to reduce the computation cost, while the image resolution and surface detection accuracy are well preserved. In fact, all these features of the proposed algorithm are not limited to endoscopic airway OCT, it could also be used to improve the image segmentation of other endoscopic OCT systems, such as intravascular OCT [35] or endoscopic gastrointestinal tract OCT [36]. The limitation of the proposed algorithm is that the detected optimal surface found by our algorithm must be continuous in the volumetric data, and thus making it difficult to process airway OCT data with bronchial bifurcations. Lastly, further improvement of the proposed algorithm is possible by incorporating edgeenhanced weights in the 3D graph, such as using better edge representations trained from a convolutional neural network [37]. To summarize, we proposed an automatic 3D segmentation method to identify and detect multiple tissue layers in endoscopic airway OCT images. This three-dimensional segmentation scheme preserves the smooth and continuous airway tissue surface by performing a 3D optimal graph search constrained by locally adaptive tissue slope. The robustness of our method against multiple artifacts and noise in the OCT images was demonstrated in the segmentation of three animal and human data sets acquired by two different OCT systems. Quantitative comparisons on results obtained from the proposed method, the 2D-based DP method, the conventional 3D GS method and the gold-standard manual-labeled method further verified the accuracy and superior performance of the new approach.