Automatic correction of the initial rotation angle error improves 3D reconstruction in endoscopic airway optical coherence tomography

: Endoscopic airway optical coherence tomography (OCT) is an advanced imaging modality capable of capturing the internal anatomy and geometry of the airway. Due to fiber-optic catheter bending and friction, the rotation speed of the endoscopic probe is usually non-uniform: at each B-scan image, the initial rotation angle of the probe is easily misaligned with that of the previous slices. During the pullback operation, this initial rotation angle error (IRAE) will be accumulated and will result in distortion and deformation of the reconstructed 3D airway structure. Previous attempts to correct this error were mainly manual corrections, which are time-consuming and suffered from observer variation. In this paper, we present a method to correct the IRAE for anatomically improved visualization of the airway. Our method derived the rotation angular difference of adjacent B-scans by measuring their contour similarity and then tracks the IRAE by formulating its continuous drift as a graph-based problem. The algorithm was tested on a simulated airway contour dataset, and also on experimental datasets acquired by two different long range endoscopic airway OCT platforms. Effective and smooth compensation of the frame-by-frame cost measured the similarity between adjacent B-scans. In both simulation and experimental acquired airway OCT datasets, the algorithm achieved accurate IRAE correction and significantly improved 3D visualization of the correct airway anatomy. The proposed algorithm is computationally fast and has the potential to be directly applied to clinical airway OCT imaging applications.


Introduction
Endoscopic optical coherence tomography (OCT) is a high-resolution imaging modality that utilizes nonionizing near-infrared laser to perform cross-sectional imaging of various internal organs such as coronary arteries [1,2], gastrointestinal tract [3,4], and airway [5,6]. Incorporating a miniature rotary fiber-optic probe with linear pullback, OCT can be used to image airway structures such as mucosa and sub-mucosa with a penetration depth of ∼2 to 3 mm. Because of its unique micro-tomographic imaging feature, has been successfully evaluated in pre-clinical [7,8] and clinical [9,10] imaging scenarios.
To facilitate lateral scanning, there are two types of rotary fiber optic systems for scanning and imaging the airway: proximal scanning system and distal scanning system. In proximal scanning system [11,12], the endoscopic probe must be equipped with a torque applying mechanism (such as a torque coil) to translate rotation from the proximal end to the imaging optics at the distal end. Because the airway is a curved tubular structure, the insertion and placement of the catheter will inevitably cause the bending of the torque coil that affects rotation transmission. Moreover, the plastic protective sheath used to isolate the probe from tissue surface will be in contact Both NURD and IRAE are problems of different dimensions caused by non-uniform rotation speed. Previous approaches to correct for NURD usually treat the IRAE and NURD correction problems as a whole, which increases the complexity of the correction problem and sacrifices algorithm robustness. Also, attempting to correct for the IRAE in airway OCT using traditional NURD correction methods faces additional challenge because most of these methods relied on analyzing the texture (or speckle pattern) of the imaged tissue. However, in airway OCT, due to the long imaging range, image texture is usually lacking [as shown in Fig. 1 (b)]. Also because of the long imaging range, the elongated optical focus degrades the lateral resolution of the image, making texture-based image analysis even more difficult. Hardware correction of the IRAE may be achieved by attaching markers on the protective sheath as in gastrointestinal balloon-catheter OCT system [17]. However, balloon-catheters are not suitable for airway imaging, and the sheath markers will inevitably block a large portion of the FOV.
To date there is no previous attempt to treat NURD and IRAE corrections as two separated problems and correct for IRAE solely by using image post-processing methods. In fact, in clinical applications, IRAE correction are mostly performed by manual adjustment [19]. This manual correction is carried out by visually comparing the orientations of adjacent B-scans frame by frame. Therefore, for a regular airway OCT dataset that contains hundreds or thousands of images, this procedure is often time-consuming. Since manual operations are subjective, the alignment results may also suffer from intra-and inter-observer variation.
Addressing the above challenge, we present an automatic algorithm to correct for the IRAE of endoscopic airway OCT images. To avoid the texture missing problem, our method takes the airway contour detected from the OCT images as input, and derives the rotation angular difference between adjacent B-scans by measuring their contour similarity. With the graph-theory-based dynamic programming (DP) algorithm, the method attempts to find an optimal path that represents IRAE drifting along the pullback direction. The proposed algorithm was tested on three different datasets including a simulated dataset, a sheep airway dataset, and an upper airway dataset of a healthy male adult. The performance of our new method was evaluated by qualitative and quantitative comparison to manual aligned correction, and the results showed that our algorithm can achieve accurate, robust, and fully automatic IRAE correction for airway OCT images.

Algorithm overview
Our algorithm is designed based on two observations: 1) airway OCT lacks image texture, and 2) IRAE drifting is usually slow, continuous and accumulated. For the former, we propose to use airway contour to measure the IRAE between adjacent B-scans instead of image texture. For the later, we use a graph-based optimal path searching algorithm to enforce continuity during the tracking of IRAE drifting. The whole flowchart of the proposed IRAE correction algorithm is shown in Fig. 2. The inputs to our algorithm are the raw OCT image stack and the extracted airway lumen contour. The contour can be either manually labelled, or detected by an automatic computer algorithm such as the ones we proposed previously [20,21], the performance of which have been validated in pre-clinical and clinical airway OCT imaging applications [22][23][24].
The algorithm is divided into three steps: pre-processing, cost matrix construction, and optimal path searching. In the pre-processing step, the airway contour is smoothed and resampled to exclude the shape difference in the airway contours between adjacent B-scans, thereby ensuring the accuracy of the subsequent IRAE correction. In the cost matrix construction step, a cost matrix is constructed through measuring the similarity of adjacent B-scans. In the optimal path searching step, the graph-theory-based dynamic programming (DP) algorithm is performed on the cost matrix to search for the optimal path representing the drifting of IRAE. Finally, the corrected OCT images are obtained by aligning the images using the estimated IRAE value.

Pre-processing
The purpose of the pre-processing step is to eliminate shape-induced contour errors and ensure reliable similarity measurement. It includes the following procedures: 2) Contour smoothing: each airway contour is smoothed with a moving average filter to remove part of the outliers.
3) Contour resampling: resample every airway contour at a selected sampling rate to mitigate the lumen shape error.
The most important procedure in the pre-processing step is contour resampling. To access the degree of IRAE between two consecutive B-scans, we can measure the similarity of the two airway contours (e.g. using the L2 norm). Ideally, the similarity between adjacent B-scans should only be originated from the inconsistent initial rotation angular position. However, the shape and geometry of the airway lumen contour is continuously varying along the pullback direction. Part of the airway may be even missing due to limited FOV. In addition, the accuracy of the detected contour worsens when the obtained images are contaminated with noise or motion blur. As shown in Fig. 3 (a-b), the contour shape is obviously different on two adjacent B-scans [yellow dashed box in Fig. 3(b)] due to detection error of the airway lumen. The difference caused by the shape of the adjacent airway contours will adversely affect the calculation of the contour similarity, resulting in a wrong IRAE correction result. The resampling operation is to address the above problem. Its principle and operation are detailed as follow: 1) For each B-scan, sort the contour points in the current B-scan according to the point difference d k+1,i : where C k,i is the Cartesian coordinate of the i-th point on the airway contour of the k-th B-scan, and m is the total number of B-scans.
2) Based on the sorted data, resample the current contour by selecting N resample contour points with less than d k,i value according to the sampling rate s: where s ∈ [0, 100%], N total is the total number of contour points.
3) Calculate the Cost between the resampled contour and the reference contour (previous B-scan): where C k,i is the Cartesian coordinate of the i-th point on the airway contour of the k-th B-scan, N resample is the number of points on the resampled contour, and I resample is the resampled contour point set. By looping through different s, a Cost-to-s diagram is obtained for each B-scan, as shown in Fig. 3 (c).  Equivalent to perform a threshold segmentation on Fig. 3(d), the Cost threshold is either user defined or automatically selected. Figure 3 (e) shows the result of the Cost segmentation of Fig. 3(d) using a threshold obtained automatically by the Ostu's method [25]. The yellow dots in Fig. 3(e) represent the specific sampling rates determined for each B-scan, and Fig. 3 (f) shows the contour resampling result of Fig. 3 (b) at the selected sampling rate. In addition, the Cost threshold can be further fine-tuned for each dataset to obtain the best result, thus ensuring the accuracy of resampling operation.
After all the B-scans are resampled, a resampled contour set is obtained, which is used to calculate the IRAE at the next step. The above resampling operation allows us to use only one parameter (the Cost threshold) to mitigate contour error of the whole dataset. An effective pre-processing provides a solid foundation for the accurate correction of initial rotation angle.

Cost matrix construction
The correction of IRAE is equivalent to rotate each polar image clockwise or counterclockwise by an estimated angle, and the counterclockwise (clockwise) rotation in polar coordinate system is equivalent to the right (left) cyclic shift in Cartesian coordinate system. For example, +n represents a right cyclic shift of n units in Cartesian coordinate system, which correspond to a counterclockwise rotation of 2πn/N total radian in polar coordinate system. Therefore, to estimate the IRAE at each B-scan, we can first shift the current contour for different units, and then compare each shifted contour to the contour in the previous B-scan. The comparison result is then store into a matrix with dimension (2n+1) × (m-1), which we refer to as the cost matrix C. In this way, each element in the cost matrix represents the similarity between a shifted contour and the reference contour. Figure 4 shows the flow diagram of cost matrix construction. The original un-processed contour set is referred to as the reference contour set. To construct the cost matrix, the airway contour that has been resampled is first cyclically shifted by -n to + n units, and then the Cost between the shifted contour and its reference contour is obtained after each shift. Notice that the value of n should be chosen larger than the maximum initial rotation angle deviation between two contours. Due to the resampling and cyclic shift operations, the element C(j,k) of the j-th row and the k-th column in the cost matrix is modified from (3) to: where C j k+1,i is the Cartesian coordinate of the i-th point on the (k+1)-th airway contour after j-th shift, which corresponded to cyclically shifting by (n-j) units. As shown in Fig. 4, the red numbers in the cost matrix C indicate the cost values obtained by (4), and a global optimal path (black curve) representing IRAE drifting is identified from left to right by an optimal path searching algorithm.

Optimal path search
Because the variation in the initial rotation angle of adjacent OCT slices is generally small, and its changing rate is slow and continuous during the pullback, the drifting of IRAE should be a smooth trajectory. Therefore, an optimal path searching algorithm that imposed continuity constraint during searching should be used for IRAE estimation on the cost matrix. In this work, we use the dynamic programming (DP) algorithm [26], which is a graph-theory-based optimal path searching technique, to meet the above continuity requirements. The concept of DP algorithm is to decompose the original problem into multiple sub-problems, and then recursively solve the original problem at O(n 2 ) complexity. Due to its fast and stable performance, DP algorithm has been used to solve various optimization problems including gaming designs [27], economic planning [28], computer vision [29,30] and etc.
The continuity constraint in the DP algorithm is imposed by connecting the neighborhood nodes with edges. In the original DP algorithm, the path elements are exactly one row apart in adjacent columns, that is, each element is only linked to three possible successors in the next column. In this case, some abrupt changes in the initial rotation angle will adversely affect the search for the optimal path. We proposed to modify the continuity constraint by changing a parameter p, which adjusts the search span between columns. As shown in Fig. 5, when the search span does not exceed the upper or lower boundary of the matrix, the search span is 2p+1. If the search touches the upper or lower boundary of the cost matrix, the search in that direction will be aborted. Figure 5 (b) illustrates the graph structure between two adjacent columns for p = 1, 2, and 3. A larger p value allows for better flexibility and adaptability, but increases the risk of miss estimation due to noise interference. Next, to achieve back-tracing, we calculate and assign a weight value to each node. In our implementation, the weight g(j,k) of the node in row j and column k is given by: where imposed the continuity constraint of node (j, k), and C(j,k) is the element value of the row j and column k on the cost matrix. Since g(j,1) is known, the optimal path is traced back starting from the node with the minimal weight in the last column. Finally, the rotation correction angle A k of the k-th OCT image is computed using the following equation: where P q is the row index corresponding to the q-th column on the optimal path.

Simulation setup
We constructed a simulation dataset with pre-defined IRAE to evaluate the IRAE correction performance of our proposed algorithm. We employed a customized circular contour with four orthogonal protrusions as the initial contour, as shown in Fig. 6 (a) and (b). The initial contour was transformed into 500 sampling points in the Cartesian coordinate system, and Gaussian white noise was randomly added with variance σ c . To simulate the IRAE drifting, the rotation angles during pullback were calculated according to the Logistic function. The maximum value and steepness of the Logistic curve were set to 50 and 0.1, respectively. Then we added Gaussian white noise with variance σ a to the Logistic function curve. The final IRAE drifting curve is shown as the green curve in Fig. 6 (c). A total of 100 contours were generated according to the IRAE curve.

Experimental endoscopic airway OCT datasets
We evaluated the feasibility and flexibility of our automatic IRAE correction on a sheep airway dataset and an adult human airway dataset acquired by two different long range endoscopic airway OCT imaging systems [31,32]. These long range OCT (LR-OCT) systems were designed to have an improved imaging range, either by using light source with long coherence length, or by solving the superposition interference signal via phase modulation. With the extended imaging range, LR-OCT is able to image the airways of large animals or human adults, which have internal diameters exceeding the typical OCT axial imaging range. This study had been approved by Southern Medical University and carried out in compliance with institutional guidelines. The characteristics and details of the subjects and animals in our study could be found in Table 1. To obtain the airway lumen contour, the sheep airway was manually segmented frame by frame. 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 deleting specific control points until the annotator is satisfied with the result. The human airway contours were automatically segmented on an endoscopic OCT dataset containing 2500 images of the upper airway. The automatic segmentation was done using the graph-theory-based algorithm presented in [20], which worked by recursively finding the shortest-path (with minimal cost) on the node-weighted graph derived from the image. In order to enlarge the testing dataset, we have also created several more datasets by downsampling the number of B-scans at equal intervals. For the human dataset, by setting the down-sampling intervals to 2, 5, 10, and 25, we obtained four more datasets with 1250, 500, 250, and 100 total frames respectively. And for the sheep data, we obtained one more dataset by down-sampling at 2-frames interval. We then performed IRAE correction using the proposed method on each dataset independently.
We further tested the Normalized Cross-Correlation (NCC) algorithm [33] for performance comparison on each dataset. The specific steps of the NCC algorithm are: 1) Convert the current B-scan from polar coordinates to Cartesian coordinates.
2) The border of the current B-scan is padded by 10 pixels to the left and right by replicating the boundary pixels.
3) Compute the NCC value between the previous B-scan and the padded B-scan according to: where f is the previous B-scan and n is its total number of pixels, and g is the padded B-scan. 4) Find the position corresponding to the maximum NCC value, which indicated the IRAE.

Manual IRAE correction
Manual alignment of the OCT images was performed to correct for the IRAE. This manual alignment was done in Amira (FEI Visualization Sciences Group, Burlington, MA) by an experienced annotator. During the alignment, two consecutive images were superimposed in a semi-transparent format within Amira to display tissue contours. Then the user manually adjusted individual image slices according to the similarities of the grey intensities until adjacent slices overlapped. In the case of misalignment, two different tissue contours were displayed that allowed the user to rotate either image until tissue overlap was achieved. For faster processing, the OCT images were resized to 512 × 512 pixels before alignment using IrfanView (Irfan Skiljan, Wiener Neustadt, Austria).

Evaluation metric
To assess the IRAE correction results, all the uncorrected and corrected airway contours were converted into circumferential images, and then imported into Mimics (Materialise, Leuven, Belgium) to perform 3D reconstruction. And the mean absolute deviation (MAD) and the Pearson Product-Moment Correlation Coefficient (PPMCC) are used as the validation metrics for quantifying the correction accuracy of two automatic algorithms, which are given by: where A is the auto-correction result, and M is the preset rotation angle in the simulation experiment or the manual correction result in the test experiment. For statistical analysis, differences between manual and computer corrected angles were determined using Bland-Altman plots [34] after IRAE correction. Figure 6 (a-c) show the IRAE correction results on the simulation dataset when σ c = 10 and σ a = 6. The overlapping contours after IRAE correction in the polar coordinate system and the Cartesian coordinate system are shown in Fig. 6 (a) and Fig. 6 (b), respectively. As can be seen in these figures, compared to the non-corrected contours, the misalignment of the contour protrusion has been significantly reduced after performing IRAE correction. Figure 6 (c) shows the reference, noisy, and estimated rotation angles at different frames. The rotation angles estimated by the proposed algorithm were closely in accordance with the reference rotation angles. Next, we investigated the effects of the noise parameters σ c and σ a on the IRAE correction performance. Figure 6 (d) and (e) show the MAD and PPMCC values for the algorithm evaluation under the influence of different σ c and σ a . As can be seen, the algorithm achieved error-free angle correction when σ c = 0. As σ c increased, the MAD increased and the PPMCC decreased, but both for only small amounts (MAD maintained low and PPMCC>0.99). The MAD and PPMCC values vary with different σ a , but its amplitude was also within a small range. This illustrated a robust anti-noise performance of the proposed algorithm.

Sheep airway results
The 3D reconstruction results of the sheep airway before and after IRAE correction are shown in Fig. 7 (a). As can be seen, in the uncorrected 3D model, the airway structure was distorted with severe twisting. Also, the NCC algorithm failed to recover the twisting caused by continuous IRAE, and the airway shape is distorted. In our corrected model, the twisting has been corrected for, and the structure of the airway was significantly improved with the proposed method. Moreover, the circumferential OCT image slices at the right panel of Fig. 7(a) are corresponding to the red, blue, and green circles on the 3D airway model at the left panel. As can be observed from these circumferential images, difference in the image orientation between the manual result and the result of the proposed algorithm was slight. Figure 7 (b) shows the correction angles with respect to the slice index obtained by manual correction and two automatic methods. It can be seen that the rotation angles of our algorithm correction and manual correction were roughly consistent, whereas the NCC curve deviated with the manual correction result. Moreover, a bias of 0.01236π was obtained from the Bland-Altman diagram shown in Fig. 7 (c). This means that the rotation angles determined by our method agreed very well with the manually labelled angles. Next, the MAD and PPMCC values were calculated to evaluate the performance of different algorithms on down-sampling datasets. The quantification results and computational cost of the two comparing methods are listed in Table 2. The algorithm correction result with lower MAD and higher PPMCC is more consistent with manual correction. It can be seen that compared with NCC, the correction results of our proposed method were closer to the manual correction results at all two sampling intervals.

Human airway results
The 3D reconstruction results of the human airway are shown in Fig. 8 (a), and the auto-manual correction comparison and Bland-Altman diagram are shown in Fig. 8 (b) and Fig. 8 (c) respectively. It can be seen from the 3D reconstruction results, there was no obvious structural difference between the 3D reconstructed airway models obtained through the proposed method and manual correction, while airway twisting is still visible in the NCC correction result. The estimated IRAE drifting curves in Fig. 8(b) also confirmed the poor performance of the NCC algorithm. Because the human airway dataset consisted of 2500 B-scans, to reduce the manual correction time, the annotator only corrected 250 selected key B-scans. This resulted in the jagged manual correction angle in Fig. 8 (b) and the line-segment-like distribution of the Bland-Altman diagram in Fig. 8 (c). A bias of 0.01502π was obtained from the Bland-Altman diagram, and the MAD and PPMCC calculated for the human contour set were 0.0270π and 0.9500 respectively, which were statistically significant. The auto-manual correction accuracy was slightly decreased compared to the sheep airway data. It may be because the human airway contours were automatically detected, making its accuracy not as high as the manually labelled sheep airway contours. Table 3 lists the quantification results and computational cost of the human dataset. The performance of the proposed method decreases slightly when the down-sampling interval increases due to the reduced correlation between adjacent contour features. However, our algorithm's results were still highly in consistent with the manual correction results across all different down-sampling rate. The above results successfully demonstrated the accuracy and robustness of the proposed algorithm on human airway OCT imaging.

Computational cost
In addition to the accuracy and robustness discussed above, our proposed algorithm is also fast: the theoretical computational cost of the DP algorithm is only within quadratic time. In our dataset, the airway contour at each B-scan of both the sheep and human airway datasets contains 500 points. With MATLAB-based codes running on an Intel Core i5 desktop computer, the average processing time for the sheep dataset with 220 images and for the human dataset with 2500 images were 0.3248 s and 3.6863 s respectively. Comparing the time consumption in Table 2 and Table 3, the proposed algorithm is much faster than the NCC algorithm. The mean processing time for each B-scan is 0.0015 s across all datasets. Moreover, because only the contour data from adjacent B-scans are used for alignment, the optimal path searching can be updated online as new contours are added, making real-time processing feasible.

Discussion
The non-uniform rotation effect in endoscopic airway OCT has long been a problem that degrades image quality and affects airway disease screening. In this work, we tried to correct for one of distortions caused by this non-uniform rotation effect, namely, the initial rotation angle error, or IRAE. In view of the continuous drifting feature of the IRAE, we presented a fully automatic correction method based on the dynamic programming algorithm. Our results on both simulated dataset and experimental airway datasets have shown that the proposed algorithm can achieve accurate IRAE correction across different airways with different diameters, and thus significantly improves 3D airway reconstruction. The contribution of our work includes the following: 1) Although IRAE and NURD are both originated from the non-uniform rotation speed of the imaging probe, we have demonstrated that IRAE correction could be solved by considering it as a separated problem independent of NURD.
2) Unlike most NURD correction methods, which usually rely on the texture (or speckle) information in the images, the proposed IRAE correction method is based on airway contour information. In this way, IRAE estimation accuracy can be preserved even when tissue texture is obscured (a common problem in long range endoscopic OCT).
3) The similarity of airway contours between adjacent B-scans allows IRAE to be measured, making IRAE drifting identifiable. We have imposed several measures to ensure robust IRAE estimation, namely the resampling operation to isolate contour error, and the adjustable continuity constraint in optimal path searching to accommodate contour variation.
4) The use of airway contours for IRAE estimation also lowers the computational complexity. Combined with our DP-based algorithm, real-time processing is achievable.
The proposed algorithm is based on pure post-processing of the OCT images and thus is not relevant to the OCT hardware setup. This means that our method may also work with other kinds of endoscopic OCT systems [35,36]. In addition, by optimizing the code and migrating it to other efficient platform such as C++, our method has high potential to be translated into real-world clinical applications.
In practice, when the probe stops pulling back or when the torque is accumulated to an extreme, the torque will try to recover. This may result in reverse IRAE drifting. However, in this torque recovery process, our method still has the ability to estimate the IRAE because the algorithm does not constrain the direction of the IRAE. This means that the proposed method will be feasible as long as the drifting is continuous and the dominant distortion is IRAE rather than NURD.
Our algorithm relies on the tracking of the contour features among consecutive slices. Therefore, it requires the continuity of the imaged organ to ensure correlations between OCT slices. Increasing the sampling density in both the circumferential and longitudinal directions [37,38] may reduce the interference of contour discontinuity, thereby improving algorithm performance. But still our algorithm could not handle disconnected contours such as airway bifurcations. Another potential limitation of the current study is that for samples with more circular shapes, especially when the catheter is well centered in the airway, the correction performance of our method may be insignificant. Moreover, although this algorithm is applicable to a range of scenarios where IRAE occurs, the contour may also be affected by NURD. Therefore, performing NURD correction [39,40] prior to IRAE may have the potential to improve IRAE correction accuracy.
Finally, our method compensates for the drifting motion of the imaging probe, but does not account for the motion of the airway tissue. It's worth mentioning that the airway motion cannot be neglected due to a limited pullback speed and limited imaging framerate, and some offsets would be expected in anatomical reality. In this case our method may not be able to represent the true and natural shape of the airway.

Conclusion
The initial rotation angle error in endoscopic airway OCT imaging affects the visualization of the correct airway structure. Its correction relies heavily on manual operation. To solve this, here an automatic IRAE correction method is proposed. Based on the slow and gradual changing feature of the IRAE, our approach proposed to search for an optimal correction path within a cost matrix that measured the similarity between adjacent B-scans. In both simulation and experimental acquired airway OCT datasets, the algorithm achieved accurate IRAE correction and significantly improved 3D visualization of the correct airway anatomy. The proposed algorithm is computationally fast and has the potential to be directly applied to clinical airway OCT imaging applications. Data availability. The data that support the findings of this study are available on request from the corresponding author.