Automatic airway wall segmentation and thickness measurement for long-range optical coherence tomography images

We present an automatic segmentation method for delineation and quantitative thickness measurement of multiple layers in endoscopic airway optical coherence tomography (OCT) images. The boundaries of the mucosa and the sub-mucosa layers were extracted using a graph-theory-based dynamic programming algorithm. The algorithm was tested with pig airway OCT images acquired with a custom built long range endoscopic OCT system. The performance of the algorithm was demonstrated by cross-validation between auto and manual segmentation experiments. Quantitative thicknesses changes in the mucosal layers are obtained automatically for smoke inhalation injury experiments.


Introduction
Imaging the sub-surface structure of the airway wall has great significance for detecting abnormalities during airway injuries.Smoke exposure and inhalation risks, including thermal, toxic, and chemical injuries, result in airway hyperemia, edema, sloughing and necrosis [1,2].Pathophysiological information of the injured airway wall, such as the thickness of the mucosal, the gathering of mucus, and the deformation of the airway lumen, could provide better diagnosis of respiratory decrease [3][4][5].As a nonionizing, non-invasive medical imaging modality, optical coherence tomography (OCT) has been used to perform high resolution, cross-sectional imaging of biological tissues.Previous studies [6][7][8][9][10][11] have demonstrated the flexibility of OCT to image the airways of different animals in vivo.Among them, the Fourier domain long-range swept source OCT (LR-SSOCT) exhibits higher sensitivity, faster imaging speed and an extended imaging range compared to conventional OCT systems.Figure 1(a) shows the schematic diagram of the LR-SSOCT.The wavelengthsweeping laser is split by a coupler into the sample and reference arms.Light comes back from both arms to another coupler, and the interference signal is then picked up by a photosensitive detector.A phase modulator (electro-optic or acousto-optic), which provides a frequency shift in the interferometer signal, enables an extended imaging range by preventing the superposition of the positive and negative terms [12].The sampling of the LR-SSOCT is performed by a miniature endoscopic OCT probe that can simultaneously rotate and be pulled back.As shown in the schematic drawing in Fig. 1(b), the probe is constructed by carefully aligning and splicing a 45-deg rod mirror, a GRIN lens, a spacer and a single-mode fiber within a metal housing.The length of the spacer (a portion of no core fiber) determines the working distance of the probe.The metal housing is firmly soldered to a stainless torque coil which translates the torque from the proximal end of the probe to the distal end of the probe.Two-dimensional (2-D) cross-sectional images are reconstructed from the interference signals for each rotational period.Figure 1(c) and 1(d) show the circumferential and Cartesian OCT images of one cross-section in the sheep airway, respectively.These sample images were acquired by our recently reported improved LR-SSOCT system [10,11] that can achieve an extended imaging range of 25 mm and a resolution of 10 μm in tissue.The system utilized a 1310-nm swept laser source with 50 kHz A-line rates.The imaging probe has an outer diameter of 1.47mm and rotates at 1500 rpm to achieve 25 frames/s (2000 A-lines per frame).Figure 1(e) displays the enlarged detail of the red box in Fig. 1(d).The structure of the airway wall, including the cartilage, the mucosa, the submucosa, and the mucus, are depicted clearly.The identification of different structures in the OCT images are of great importance for quantitative evaluation of airway injury.However, although we have already obtained high resolution airway wall OCT images, the classification of the airway wall structures are based on manual labeling of the boundaries [13], which is time-consuming and subject to observer errors.A clustering algorithm is used to segment the airway structure from the background [14], but the airway layers are not detected.A morphological and thresholding method is used to identify these layers [15]; however, their method is not robust enough to handle speckle noise and missing features in quality degraded images.No other method has been reported on the auto-segmentation of airway wall structures on OCT images other than these two.
In order to robustly segment different airway wall structures and provide quantitative information of these structures automatically, we present a graph theory based segmentation algorithm using Cartesian airway OCT images as input.This algorithm works in two steps, the pre-processing step, which localizes the potential airway wall regions, and the edge delineation step, which detects the exact locations of three edges specifying the airway lumen, the mucosa/submucosa boundary, and the submucosa/cartilage boundary respectively.The quantification of the average layer thicknesses can be obtained afterward.This algorithm was tested in the sheep and pig airway images acquired by our LR-SSOCT system [11].The results show that our algorithm can achieve accurate, robust, and fully automatic delineation of multiple structures in airway OCT images.

Method
The whole algorithm workflow is shown in Fig. 2. To start with, the purpose of the preprocessing step is to distinguish the entire airway structure from the background.It is mainly comprised of denoising and binarization of the original input images, and then identifying the prominent OCT airway regions.After the pre-processing step, an edge detection process where the actual airway lumen, the edge between the mucosa and submucosa, and the boundary of the airway structure are localized using a dynamic programming based approach.Finally, the thickness measurement process is carried out on the segmented airway regions where the laser beam has an orthogonal incident angle.

Pre-processing
To achieve automatic airway wall segmentation, the first step should be to identify the entire airway structure in the given OCT image.This means that the airway region should be distinguished from the background for further processing.However, the OCT images are generally distorted by multiple defects.These defects include: speckle noise inherited from laser imaging [16,17]; mirror image/objects induced by Fourier transformation [18]; and ghost objects produced by internal interference of the optics.Also, since the imaging probe is protected from the tissue by a plastic sheath, this sheath needs to be distinguished from the airway structure.As can be seen in the sample input image in Fig. 3(c), although the airway structure seems prominent by human judgment, it is still challenging for the computer to precisely identify the tissue since the aforementioned noise is also obvious.
To address these problems, we applied a series of low-level operations on the OCT input images to produce classified airway regions.These operations include: 1) Speckle noise suppression: median filtering is applied.
2) Horizontal noise reduction: the mean intensity of each horizontal line is subtracted from each pixel in the same line to suppress the horizontal artiest such as saturated lines and the sheath.
3) Binarization: the image is transformed to a black and white (BW) image by a user defined threshold value [Fig.3(b)].
5) Dilation and bridging: the candidate regions are further dilated and bridged.Some of the previously separated structures will be connected again.An example is given in Fig. 3(d).Moreover, when imaging a long airway segment, the bifurcations of the airway will sometimes be seen on the OCT frames.In the Cartesian coordinates, the airway structure will break into two disconnected regions.To deal with this, in the area filtering step of our algorithm, only the binarized regions with an area that is larger than our minimum area threshold will be preserved.The edge detection procedure will be carried out in these regions separately afterward.
The final segmented BW regions are considered to contain reliable airway structures.More specifically, the upper boundary of these regions could be considered as the coarse measurements of the airway lumen (the initial lumen).This operation is robust since it eliminates the defects in the images, such as the "stains," and can handle harsh noise caused by hardware deficiencies.

Edge detection based on dynamic programming
The successful localization of the airway regions provides solid foundation for the edge detection procedure where the precise boundary locations are indicated.The dynamic programming (DP) algorithm is the cored method introduced in this step.Developed from the graph theory, the major concept of this algorithm was first invented by Richard Bellman in the 1940s [19].Since then, this algorithm has found widespread applications including the shortest path calculation in computer gaming, the optimal consumption and saving in economic planning, and edge detections in computer vision.Prior to our work, the DP algorithm has been successfully applied to the segmentation of OCT images such as the delineation of the multiple retinal layers [20][21][22][23] and corneal layer boundaries [24].The major feature of this algorithm when applied to edge detection is that it can preserved the continuity of the boundary, which means that it is less affected by outliers.Also, compared to other #247966 graph-based edge detection methods [25], the DP algorithm is an efficient graph solving method that has a quadratic time complexity.
Generally, the concept of the DP algorithm is to break down the original problem into multiple small problems and then solve the original problem "recursively."Here, the principle of the edge detection by DP algorithm is reviewed briefly.As illustrated in an example in Fig. 4, the input image is a gradient image of the original image with the edge features protruding from the background; each pixel in the image is considered a node in a graph and thus the edge detection process is transformed to a path finding process.Generally, two steps are taken to find the optimal path: graph construction and recursive solution finding.The graph construction step is to assign cost values to each node in the graph; in our implementation, the cost for node (i, j) is given by: where λ is a scale factor; i, j are the transverse direction and vertical (depth) direction, respectively; and Calculating the cost for each node is actually measuring its similarities against three of the neighboring nodes at the right side.Thus, the graph construction process is starts from the left side of the graph and ends at the right side.While getting the cost for each node, anther parameter is also recorded for each node: which logs the possible solutions for getting the shortest path.After constructing the graph, the solution finding stage is carried out by first finding the minimum cost at the last column, which in fact represents the total minimum cost for the shortest path in the constructed graph and then locates the shortest path by retrieving the neighboring node locations recorded by the path (i, j) parameters one by one.As illustrated in Fig. 4 (b), the total minimum cost is given by the fourth node in the last column, and then the solution is resolved by interpreting the path values starting from this node and ends in the first node at the first column.In our application, the gradient image is first generated to highlight the edges in the original image.Unlike retinal OCT images, the airway structure in Cartesian OCT images is curvy.Thus, in order to get better represented edges, the gradient image G is given by the L 2 norm of the horizontal gradient g x and vertical gradient g y : Next, the detection of multiple layer structures in the airway OCT images follows three steps: 1) Graph construction: identify the region-of-interest (ROI) for each edge.To extract an unknown edge, a known reference edge will be used.Given this reference edge, two offsets are used to define the ROI for the current edge.These offsets are the positive and negative y-offsets with respect to the reference edge.They define how many pixels above and below the edge should be taken as the ROI/graph.In this way, a region of constant height could be defined, and the graph is constructed from the gradient image of the same region after flattening.
2) Edge localization: using the dynamic programming algorithm, the precise position for each edge is extracted by solving the shortest path problem in the constructed graph.
3) Edge refinement: the edge is further smoothed by averaging neighboring edge pixels.
The precise lumen is extracted as the first step.The initial lumen extracted in the preprocessing step was used as the reference edge, and a graph [Fig.5(a)] is constructed using the two offsets defined for the lumen edge.After solving the precise lumen (the green edge in Fig. 5), another graph [Fig.5(e)] was constructed using the precise lumen as the reference and two different offsets to find the outer boundary of the airway (the yellow edge in Fig. 5).After these two edges, a closed region wrapped by the lumen and the outer boundary was defined as the final graph [Fig.5(c)] to extract the boundary between the mucosa and submucosa layers, namely the middle edge, and thus the two other edges will be avoided.Figure 5 shows the intermediate results for one sample image.Figures 5(a), 5(c) and 5(e) show the flattened ROI (graph) for the lumen, the middle edge and the outer boundary, respectively.Moreover, around the bifurcation segment in the airway, the graphs for each disconnected regions are constructed separately according to the disconnected initial lumen, and the DP algorithm is carried out to subsequently find the edges for each region.
Figures 5(b), 5(d) and 5(f) show the extracted minimum path solved by the dynamic programming algorithm.Figure 5(g) shows the original airway image overlaid with the detected edges.It can be seen that the mucosa and submucosa layers are detected accurately.The smoothness of these edges is successfully preserved against multiple distractions, such as speckle noise and cartilages.
Finally, to accurately quantify the airway layer thickness, only the regions that are penetrated by the laser beam with an orthogonal direction should be taken into consideration.Therefore, after the delineation of the whole airway structure, these regions were selected based on a given flatness threshold on the Cartesian images.The average thickness of different layers could be obtained easily afterward.

OCT system setup and animal preparation
To verify the flexibility of our algorithm, the images acquired by the LR-SSOCT system reported by our group previously [10,11] were used as the test set.The LR-SSOCT system utilized a customized imaging probe that rotates at 1500 rpm with a pullback speed of 12.5mm/s.In order to achieve long range imaging, the probe was designed to have an extended working distance of 20 mm and the axial resolution was 10 μm in tissue.Imaging of airways with a maximum diameter of 50 mm and 20 cm length can be achieved with this system.
Specifically, airway OCT image data sets of two different animals acquired by this LR-SSOCT system were used: a sheep airway data set and a pig airway data set.These images were taken on convenience samples from an ongoing study which involves a clinically relevant model of lung failure due to inhalation of wood bark smoke and cutaneous burn.For #247966 both the sheep and the pig samples, baseline airway OCT images were acquired after induction of anesthesia to the animals.Additionally, for monitoring the conditions and progression of smoke inhalation injury, OCT images were acquired after smoke inhalation injury was induced [26].
Figures 6 and 7 illustrate the edge detection results in different frames in baseline data sets of the sheep and pig samples, respectively.As can be seen, even when some of the airway wall was not in the OCT imaging zone with the highest sensitivity, the airway wall was detected accurately.Also, the algorithm could effectively deal with the bifurcation of the airway; the layers in the disconnected areas were extracted successfully.

Auto-manual segmentation comparison
Manual segmentation was performed to evaluate the performance of the proposed autosegmentation method on the baseline OCT airway images.Specifically, the manual segmentation was done by clicking 30-50 points on an edge and then spline-fitting these points to determine the actual edge.The spline-fitting curve could be modified by adding or removing specific control points until the annotator was satisfied with the result.In total, 100 frames of healthy sheep airway OCT images picked randomly from data sets containing 400 images, and 50 frames of healthy pig airway OCT images picked from data sets with 200 images were manually annotated.Among the results, all three edges of interest were detected.
To quantify the measurement accuracy of the proposed segmentation algorithm against manual segmentation, two validation metrics are used: the root mean square error (RMSE) and the mean absolute deviation (MAD), which are given by: 2 1 1 1 ( ) where A is the auto-segmentation result, M is the manual-segmentation result, and n is the total number of pixels.We compared the auto-manual segmentation accuracy in the airway regions with clear layer structures, and the results are listed in Table 1.Compared to the penetration depth of the LR-SSOCT system, which is about 2 mm, the deviation of the auto and manual segmentation results is less than 1% of the total tissue depth.The error rates of the sheep and the pig airway are similar to each other, demonstrating the proposed algorithm's robustness across different animals with different airway diameter.The localization of the airway lumen and the mucosa/submucosa edge are more accurate than that of the submucosa/cartilage edge; this is believed to be caused by the signal degradation in deeper tissue.

Airway thickness changes during smoke inhalation measured by the proposed algorithm
The diagnosis of inhalation injury is a primary unresolved problem in modern burn care.Using the purposed segmentation algorithm, comparable quantitative measurements of the airway layers postinjury could be obtained automatically.For the pig data set, the OCT images acquired at baseline and 1-hr post smoke injury were analyzed using our algorithm.
For each time point, 10 images were used for layer thickness measurement.These images were all acquired at the same landmark, which was the right mainstem bronchus near the proximal secondary bronchus branch.As stated previously, only the orthogonally laser penetrating areas of the airway were selected for thickness measurement.
The detection results for the pig airway images in the same location at baseline and post smoke are shown in Fig. 8.The magnified regions display the automatically selected regions for thickness quantification.The circumferential images (transformed from the Cartesian images below) in the right column of Fig. 8 provide better views of the airway substructures and the segmentation results.To verify the thickness measured by the automatic algorithm, manual results were obtained on all the images with the same protocol described in Section 3.2.Here we use three quality metrics for comparison: Average Thickness (AVGT), Root-Mean-Square Thickness Difference (RMSTD), and Dice's similarity Coefficient (DSC): where A and M are the auto and manual segmentation edges at pixel i, with the superscript t and b indicating the top and bottom edges, respectively; X and Y are the auto and manual segmented layer regions; and n is the total pixel number.Table 2 listed the comparison results for the pig data set.Even though the thickness of both the mucosa and submucosa layer increased post smoke, the RMSTD between automanual segmentation remained at a relatively low value.The segmentation accuracy for the submucosa layer suffered from a small decrease at the post smoke time point, the DSC for the auto-selected region dropped 9.6% compared to the baseline value.The sheep data set contains OCT images acquired at baseline, 1-hr, 24-hrs, and 48-hrs postinjury.These images at different time points were analyzed using our algorithm with the same set of parameters.Figure 9 shows the sample airway images and the edge detection results in the same location at different time points.The magnified regions of the Cartesian images and the circumferential images in Fig. 9 show better details of the airway substructures, the segmentation results and the automatic selected regions for thickness quantification.It can be seen that all the three edges were determined accurately using the proposed algorithm, even under obvious thickness changes across different time points.
Figure 10 illustrates the thickness changes of the mucosa and submucosa layers measured by both automatic algorithm and manual operation across these four time points.The thickness of these layers both increase as the time postinjury increases within this 48-hrs period, and the thickness of the submucosa layer grew more dramatically than the mucosa layer.Also notice that although the thicknesses of both layers changed at different time points, the values measured by the auto and the manual method were very close.The maximum deviations between the auto and manual results were less than 0.022mm and 0.21mm for the mucosa and submucosa layers, respectively.More comparison results are listed in Table 3.The DSC for all the tested images were at high values indicating that the regions segmented by the automatic algorithm were identical to those segmented manually.The RMSTD between auto and manual segmentation for layer 1 (the mucosa layer) remains below 0.06 mm across all four time points while for layer 2, the RMSTD value increases from 0.036 mm in baseline to 0.328 mm in 48 hours.This is due to the limitation of the penetration depth of the OCT system: in the postinjury case, the thickness of the airway grows almost two-fold in 48 hours compared to the baseline.The OCT signal was attenuated in deeper tissue which led to blurry edges between the submucosa and the cartilages.

Discussion
We presented, for the first time, a graph-theory-based automatic segmentation method for airway OCT image segmentation and demonstrated that this method could be implemented to successfully solve layer segmentation problems in airway OCT images.The robustness of our approach is ensured in two ways: the pre-processing is simple yet stable enough for regular OCT noise, and the dynamic programming method ensures precise extraction while preserving the smoothness of the edges.Also, our method is fast: the computational cost for the DP algorithm is within quadratic time.The average processing time for a single OCT image at the size of 2000*2000 pixels is 4.176 seconds with MATLAB based codes running on an Intel Core i7 workstation.This means that to process a full data set containing 400 images will take less than half an hour.The algorithm was based on pure post-processing of the OCT images, thus it is not relevant to the OCT hardware setup.Thus our method can be integrated to real world clinical applications directly, and the diagnostic efficiency can be greatly improved.More importantly, the automatic processing also gives quantitative measurement of pathophysiological parameters of the airway which could benefit the accurate diagnosis of airway abnormalities by given comparable standards.Compared to manual segmentation, these parameters remove any observer variations and, therefore, are more reliable.
The proposed algorithm has a number of limitations.In some cases, the outer boundary of the airway could not be successfully delineated due to the limited penetration depth and the non-uniform axial sensitivity distribution of OCT systems.In addition, since the airway structure is not completely circular and the probe is not always in the center, when the OCT beam is near tangent to the tissue surface, the OCT signal is degraded.Improving the OCT hardware performance could minimize these problems.Finally, to improve the processing speed of the proposed algorithm, the major barrier is the computational cost for solving the shortest path problem in the constructed graph.This process can be accelerated by implementing the algorithm on parallel processors such as the Graphical Processing Units (GPUs).

Conclusion
A fully automatic airway wall structure segmentation method for endoscopic optical coherence tomography images is presented in this paper.Based on the graph theory and the dynamic programming algorithm, the proposed approach could be used to delineate the boundaries among the mucosa, the submucosa and the cartilages in upper airway wall OCT images.We demonstrated the robustness of our algorithm against multiple defects and noise in the raw input image in manual-auto cross validation experiments.Both sheep and pig airway OCT images acquired during smoke inhalation were evaluated by our algorithm, and accurate quantitative thickness changes in two different layers were obtained.

Fig. 1 .
Fig. 1.(a) The LR-SSOCT setup; (b) the structure of the imaging probe for airway imaging; (c) a sample image in circumferential coordinate system; (d) Cartesian coordinate image of (c); (e) the enlarged detail of the red box in (d).

Fig. 2 .
Fig. 2. The workflow of the proposed automatic detection algorithm.

Fig. 3 .
Fig. 3.The pre-processing step: (a) original input image where the sheath, speckle noise, and mirror objects are seen; (b) the binaried black and white image; (c) the first morphological operation: picking out major components; (d) the second morphological operation: bridging and enlarging the potential airway wall regions.

Fig. 4 .
Fig. 4. DP principle: (a) a given image where the number is the pixel intensity, and the blue line is the edge.(b) The path and cost calculated according to Eq. (1).The minimum cost is given by the 4th node in the last column; the shortest path is traced back starting from this node (green arrows).

Fig. 5 .
Fig. 5.The process of finding multiple edges: (a), (c), and (e) are a flattened graph constructed from gradient images of the airway lumen, the mucosa and the submucosa layer; (b), (d), and (f) are the corresponding shortest paths obtained by the DP algorithm of (a), (c) and (e), respectively; (g) is the original OCT image overlaid with the edges detected.

Fig. 6 .
Fig. 6.Illustration of the edge detection results in different OCT airway images taken from the pig data set.From left to right: the original images, the Cartesian coordinate images with edge detected, the circumferential images.

Fig. 8 .
Fig. 8. Thickness measurement comparison between baseline and post-smoke images of the pig data set.From left to right: original images, edge extracted images (thickness measured regions enlarged), and circumferential images.

Fig. 9 .
Fig. 9. Thickness measurement comparison between baseline and post-smoke images of the sheep data set.From left to right: original images, edge extracted images (thickness measured regions enlarged), and circumferential images.

Fig. 10 .
Fig. 10.Thickness of the (a) mucosa and (b) submucosa changes according to the time postsmoke.