Open Access
Add to BibSonomy
Abstract
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 initial angle difference was achieved. Our method has real-time capability and thus has the potential to improve clinical imaging efficiency.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. 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 with the probe, and friction will be created. This also affects the smooth rotation of the probe. Distal scanning system [12–14] employing miniature rotation device (such as a micro-motor) on the probe head can reduce mechanical friction. However, the imaging field of view (FOV) in distal scanning systems is usually incomplete due to electric wire blockage. It is also difficult to maintain uniform rotation at higher speed because of the difficulty in fabricating perfectly balanced micro-motors [15].
Whether in proximal or distal scanning system, the above hardware limitations make it difficult to ensure a constant rotation speed during endoscopic imaging, causing image distortion often known as the Non-Uniform Rotation Distortion (NURD) [16]. NURD leads to stretched or compressed B-scan images because of missing or repeated A-lines. Various methods have been proposed to correct NURD in previous researches. For example, it can be corrected by non-rigid registration using fiducial markers located on the catheter [15], analyzing speckle statistics between A-lines from intensity-based dynamic light scattering [17], or maximizing the similarity between A-lines of the aligned regions in adjacent B-scans [18]. These methods successfully improved image fidelity.
Except for NURD, the non-uniform rotation of the probe also causes difference in the beginning and ending of the images. As shown in the schematic drawing in Fig. 1 (a), when the probe scans the airway to form each B-scan image, the initial rotation angle is slightly different. Figure 1 (b) shows three consecutive OCT images in a polar image stack that suffered from this initial angular position inconsistency. To correct this kind of error, frame k+1 and frame k+2 should be aligned with frame k by rotating counterclockwise. If the uncorrected OCT images are used for 3D reconstruction, the structure of the reconstructed airway model will be distorted [Fig. 1 (c)]. In this paper, we refer to this kind of error as the Initial Rotation Angle Error (IRAE).
Fig. 1. (a) Schematic diagram of the IRAE produced during OCT pullback. (b) Three adjacent images presenting IRAE: let frame k be the reference frame, frame k+1 and frame k+2 present progressive IRAE. (c) An IRAE-uncorrected 3D model of sheep airway with severe structural distortion.
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.
2. Method
2.1 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–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.
2.2 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:
- 1) Coordinate conversion: the detected airway contours are converted from polar to Cartesian coordinates.
- 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 dk+1,i: where Ck,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.
- 3) Calculate the Cost between the resampled contour and the reference contour (previous B-scan):$$Cost = \frac{1}{{{N_{resample}}}}\sqrt {\sum\limits_{i \in {{\mathbf I}_{{\mathbf resample}}}}^{} {{{({C_{k + 1,i}^{} - {C_{k,i}}} )}^2}} } ,k \in \{ 1,2,\ldots ,m - 1\}, $$where Ck,i is the Cartesian coordinate of the i-th point on the airway contour of the k-th B-scan, Nresample is the number of points on the resampled contour, and Iresample 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).
Fig. 3. Image pre-processing to exclude the shape difference in adjacent airway contours. (a-b) Coordinate conversion and contour smoothing. (c) The Cost versus different sampling rate of one B-scan. (d) Cost-to-s diagram of the whole dataset. (e) A Cost threshold can be determined by auto-segmentation of (d). (f) Contour resampling result: the original contour in (b) is resampled to exclude miss-detected edge points.
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.
2.3 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/Ntotal 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:
Fig. 4. Flow diagram of the cost matrix construction process: 2n+1 represents the tolerance range of angular deviation, and m represents the total number of OCT images.
2.4 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(n2) 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.
Fig. 5. (a) The continuity constraint during optimal path searching is enforced by the edges between each node in the current column and the nodes in its previous column, and the search span of the nodes are 2p+1. (b) Graph structures when the continuity constraint parameter p = 1, 2, and 3.
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:
3. Experimental setup
3.1 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.
Fig. 6. Simulation results. (a) All contours overlaid in the polar coordinate system before and after IRAE correction: the black contour represents the noiseless initial contour. (b) All contours overlaid in the Cartesian coordinate system before and after IRAE correction. (c) The initial rotation angles at different B-scans. (d) MAD obtained with different σc under different σa. (e) PPMCC obtained with different σc under different σa.
3.2 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.
Table 1. Testing datasets used in the experiments
In order to enlarge the testing dataset, we have also created several more datasets by down-sampling 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.
3.3 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).
3.4 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:
4. Results
4.1 Simulation results
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.
4.2 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.
Fig. 7. IRAE correction results of the sheep airway OCT images. (a) 3D reconstruction results and example circumferential images in the cases of uncorrected, manual correction, proposed method and NCC registration. (b) The correction angle at different B-scans obtained by manual correction, proposed algorithm and NCC registration: the red, blue and green dashed lines correspond to the B-scans shown in (a). (c) Bland-Altman plot between the proposed method and manual correction: Ak is the k-th B-scan rotation angle estimated by the proposed method, and Mk is the k-th B-scan rotation angle estimated by manual correction.
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.
Table 2. IRAE correction performance of the sheep dataset at different down-sampling intervals
4.3 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.
Fig. 8. IRAE correction results of the adult human airway OCT images. (a) 3D reconstruction results and example circumferential images in the cases of uncorrected, manual correction, proposed method and NCC registration. (b) The correction angle at different B-scan obtained by manual correction, proposed algorithm and NCC registration: the red, blue and green dashed lines correspond to the B-scans shown in (a). (c) Bland-Altman plot between the proposed method and manual correction: Ak is the k-th B-scan rotation angle estimated by the proposed method, and Mk is the k-th B-scan rotation angle estimated by manual correction.
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.
Table 3. IRAE correction performance of the human dataset at different down-sampling intervals
4.4 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.
5. 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.
6. 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.
Funding
National Natural Science Foundation of China (31700857); Key-Area Research and Development Program of Guangdong Province (2018B030333001); Pearl River Talented Young Scholar Program of Guangdong Province (2017GC010282); Basic and Applied Basic Research Foundation of Guangdong Province (2021A1515012542); Guangzhou S&T Program (201804010375).
Disclosures
The authors declare no conflicts of interest.
Data availability
The data that support the findings of this study are available on request from the corresponding author.
References
1. I. K. Jang, B. E. Bouma, D. H. Kang, S. J. Park, S. W. Park, K. B. Seung, K. B. Choi, M. Shishkov, K. Schlendorf, and E. Pomerantsev, “Visualization of coronary atherosclerotic plaques in patients using optical coherence tomography: comparison with intravascular ultrasound,” J. Am. Coll. Cardiol. 39(4), 604–609 (2002). [CrossRef]
2. T. Kubo, T. Akasaka, J. Shite, T. Suzuki, S. Uemura, B. Yu, K. Kozuma, H. Kitabata, T. Shinke, and M. Habara, “OCT compared with IVUS in a coronary lesion assessment: the OPUS-CLASS study,” JACC: Cardiovascular Imaging 6(10), 1095–1104 (2013). [CrossRef]
3. K. Kobayashi, J. A. Izatt, M. D. Kulkarni, J. Willis, and M. V. Sivak, “High-resolution cross-sectional imaging of the gastrointestinal tract using optical coherence tomography: preliminary results,” Gastrointestinal Endoscopy 47(6), 515–523 (1998). [CrossRef]
4. Brand, Poneros, and Bouma, “Optical coherence tomography in the gastrointestinal tract,” Endoscopy (2000).
5. N. Hanna, D. Saltzman, D. Mukai, Z. Chen, S. Sasse, J. Milliken, S. Guo, W. Jung, H. Colt, and M. Brenner, “Two-dimensional and 3-dimensional optical coherence tomographic imaging of the airway, lung, and pleura,” The Journal of Thoracic and Cardiovascular Surgery 129(3), 615–622 (2005). [CrossRef]
6. H. O. Coxson, P. R. Eastwood, J. P. Williamson, and D. D. Sin, “Phenotyping airway disease with optical coherence tomography,” Respirology 16(1), 34–43 (2011). [CrossRef]
7. S. W. Lee, A. E. Heidary, D. Yoon, D. Mukai, T. Ramalingam, S. Mahon, J. Yin, J. Jing, G. Liu, and Z. Chen, “Quantification of airway thickness changes in smoke-inhalation injury using in-vivo 3-D endoscopic frequency-domain optical coherence tomography,” Biomed. Opt. Express 170, 295–302 (2011).
8. Y. Miao, J. C. Jing, D. Vineet, S. B. Mahon, B. Matthew, L. A. Veress, C. W. White, and Z. Chen, “Automated 3D segmentation of methyl isocyanate-exposed rat trachea using an ultra-thin, fully fiber optic optical coherence endoscopic probe,” Sci. Rep.8(1), 8713 (2018). [CrossRef]
9. R. A. Mclaughlin, J. P. Williamson, M. J. Phillips, J. J. Armstrong, and D. D. Sampson, “Applying anatomical optical coherence tomography to quantitative 3D imaging of the lower airway,” Opt. Express 16(22), 17521–17529 (2008). [CrossRef]
10. J. J. Armstrong, M. S. Leigh, D. D. Sampson, J. H. Walsh, D. R. Hillman, and P. R. Eastwood, “Quantitative upper airway imaging with anatomic optical coherence tomography,” Am J Respir Crit Care Med 15, 173 (2006). [CrossRef]
11. V. Volgger, G. K. Sharma, J. C. Jing, Y. S. A. Peaks, A. C. Loy, F. Lazarow, A. Wang, Y. Qu, E. Su, and Z. Chen, “Long-range Fourier domain optical coherence tomography of the pediatric subglottis,” International Journal of Pediatric Otorhinolaryngology 79(2), 119–126 (2015). [CrossRef]
12. D. C. Adams, Y. Wang, L. P. Hariri, and M. J. Suter, “Advances in endoscopic optical coherence tomography catheter designs,” IEEE J. Sel. Top. Quantum Electron. 22, 1 (2015). [CrossRef]
13. T. H. Tsai, B. Potsaid, Y. K. Tao, V. Jayaraman, J. Jiang, P. J. S. Heim, M. F. Kraus, C. Zhou, J. Hornegger, and H. Mashimo, “Ultrahigh speed endoscopic optical coherence tomography using micromotor imaging catheter and VCSEL technology,” Biomed. Opt. Express 4(7), 1119–1132 (2013). [CrossRef]
14. T. Wang, G. Van Soest, d. S. Van, and F. W. Antonius, “A micromotor catheter for intravascular optical coherence tomography,” Engineering 1(1), 015–017 (2015). [CrossRef]
15. O. O. Ahsen, H. C. Lee, M. G. Giacomelli, Z. Wang, and J. G. Fujimoto, “Correction of rotational distortion for catheter-based en face OCT and OCT angiography,” Opt. Lett. 39(20), 5973 (2014). [CrossRef]
16. Yoshiaki, Kawase, Yoriyasu, Suzuki, Fumiaki, Ikeno, and Ryuichi, “Comparison of nonuniform rotational distortion between mechanical IVUS and OCT using a phantom model,” Ultrasound in Medicine & Biology 33(1), 67–73 (2007). [CrossRef]
17. N. Uribe-Patarroyo and B. E. Bouma, “Rotational distortion correction in endoscopic optical coherence tomography based on speckle decorrelation,” Opt. Lett. 40(23), 5518–5521 (2015). [CrossRef]
18. G. V. Soest, J. G. Bosch, and A. F. W. V. D. Steen, “Azimuthal registration of image sequences affected by nonuniform rotation distortion,” IEEE Trans. Inform. Technol. Biomed. 12(3), 348–355 (2008). [CrossRef]
19. D. Protsenko, Z. Chen, and B. J. Wong, “Constructing 3D models of the pediatric upper airway from long range optical coherence tomography images,” Proceedings of SPIE - The International Society for Optical Engineering89262L, 89262L–89212 (2014).
20. L. Qi, S. Huang, A. E. Heidari, C. Dai, J. Zhu, X. Zhang, and Z. Chen, “Automatic airway wall segmentation and thickness measurement for long-range optical coherence tomography images,” Opt. Express 23(26), 33992 (2015). [CrossRef]
21. L. Qi, K. Zheng, X. Li, Q. Feng, and W. Chen, “Automatic three-dimensional segmentation of endoscopic airway OCT images,” Biomed. Opt. Express 10(2), 642 (2019). [CrossRef]
22. Yusi, Miao, Matthew, Brenner, Zhongping, and Chen, “Endoscopic optical coherence tomography for assessing inhalation airway injury: a technical review,” Otolaryngology (Sunnyvale, Calif.) 9 (2019).
23. K. M. Kozlowski, G. K. Sharma, J. J. Chen, L. Qi, and J. F. Wong, “Dynamic programming and automated segmentation of optical coherence tomography images of the neonatal subglottis: enabling efficient diagnostics to manage subglottic stenosis,” J. Biomed. Opt. 24(09), 1 (2019). [CrossRef]
24. Y. Miao, J. H. Choi, L. D. Chou, V. Desai, and Z. Chen, “Automatic proximal airway volume segmentation using optical coherence tomography for assessment of inhalation injury,” J Trauma Acute Care Surg. 87(1S), S132–S137 (2019). [CrossRef]
25. N. Ostu, O. Nobuyuki, and N. Otsu, “A threshold selection method from gray- level histogram IEEE transactions on systems,” IEEE Trans. Syst., Man, Cybern. 9(1), 62–66 (1979). [CrossRef]
26. R. Bellman, On the Theory of Dynamic Programming (INFORMS, 1956).
27. J. Freeman, “Dynamic programming and the evaluation of gaming designs,” OR Insight 23(2), 124–135 (2010). [CrossRef]
28. Domenico, Campisi, Massimo, GastaldiAgostino, La, and Bella, “Optimal growth and planning in a multi-regional economy: A computer program and application to the Italian case,” Computational Economics (1993).
29. B. Lee, J. Y. Yan, and T. Zhuang, “Dynamic-programming-based algorithm for optimal edge detection in ultrasound images,” Proceedings of SPIE - The International Society for Optical Engineering4549, 0193 (2001).
30. N. Merlet and J. Zerubia, “New prospects in line detection by dynamic programming,” IEEE Trans. Pattern Anal. Machine Intell 18(4), 426–431 (1996). [CrossRef]
31. L. Chou, A. Batchinsky, S. Belenkiy, J. Jing, T. Ramalingam, M. Brenner, and Z. Chen, “In vivo detection of inhalation injury in large airway using three-dimensional long-range swept-source optical coherence tomography,” J. Biomed. Opt. 19(3), 036018 (2014). [CrossRef]
32. J. Jing, L.-d. Chou, E. Su, B. Wong, and Z. Chen, “Anatomically correct visualization of the human upper airway using a high-speed long range optical coherence tomography system with an integrated positioning sensor,” Sci. Rep.6(1), 39443 (2016). [CrossRef]
33. J. N. Sarvaiya, S. Patnaik, and S. R. Bombaywala, “Image registration by template matching using normalized cross-correlation,” 2009 International Conference on Advances in Computing, Control, and Telecommunication Technologies, 819-822 (2009).
34. J. M. Bland and D. G. Altman, “Statistical methods for assessing agreement between two methods of clinical measurement,” Lancet1, 307 (1986).
35. M. L. Dufour, C. E. Bisaillon, G. Lamouche, S. Vergnole, M. Hewko, F. D’Amours, C. Padioleau, and M. Sowa, “Tools for experimental characterization of the non-uniform rotational distortion in intravascular OCT probes,” Proceedings of Spie the International Society for Optical Engineering7883, (2011).
36. C. Sun, F. Nolte, K. H. Y. Cheng, B. Vuong, K. K. C. Lee, B. A. Standish, B. Courtney, T. R. Marotta, A. Mariampillai, and V. X. D. Yang, “In vivo feasibility of endovascular Doppler optical coherence tomography,” Biomed. Opt. Express 3(10), 2600–2610 (2012). [CrossRef]
37. T. Wang, W. Wieser, G. Springeling, R. Beurskens, C. T. Lancee, T. Pfeiffer, A. F. W. V. D. Steen, R. Huber, and G. V. Soest, “Intravascular optical coherence tomography imaging at 3200 frames per second,” Opt. Lett. 38(10), 1715–1717 (2013). [CrossRef]
38. T. Wang, T. Pfeiffer, J. Daemen, F. Mastik, and G. V. Soest, “Simultaneous morphological and flow imaging enabled by megahertz intravascular Doppler optical coherence tomography,” IEEE Transactions on Medical Imaging 39, 1 (2019).
39. T. Nguyen, O. Ahsen, K. Liang, J. Zhang, H. Mashimo, and J. Fujimoto, “Correction of circumferential and longitudinal motion distortion in high-speed catheter/endoscope-based optical coherence tomography,” Biomed. Opt. Express 12(1), 226–246 (2021). [CrossRef]
40. J. Mavadia-Shukla, J. Zhang, K. Li, and X. Li, “Stick-slip nonuniform rotation distortion correction in distal scanning optical coherence tomography catheters,” J. Innov. Opt. Health Sci. 13(06), 2050030 (2020). [CrossRef]
