2.2.1.

Preprocessing algorithm

Before the thresholding can be applied, we first use a $10\times 10$ Gaussian filter with a standard deviation of 4 for every B-scan image ($I_0$) to smooth the speckle pattern. Then, we remove the saturation artifacts introduced by light reflected from a highly specular surface that are over the dynamic range of the data-acquisition system. We simplified the saturation removal method based on the minimum variance mean-line subtraction described in Refs. 23 and 24, to reduce the computation cost. Specifically, $I$, the intensity value of the saturation-free image, is assigned as

Next, to facilitate connected component analysis, we convert every grayscale B-scan into a binary image, based on a threshold value. More precisely, we applied the adaptive thresholding method²⁵ to calculate the thresholding value, and this ensures that the sample surface constitutes the largest area of bright components in the image after thresholding. Figure 3(a) shows the ex vivo bovine B-scan image. Figure 3(b) shows the corneal B-scan image after this thresholding was applied and clearly shows the upper and lower boundaries of the cornea sample.

Fig. 3

OCT images showing the result of the preprocessing steps, using a bovine corneal sample. (a) Original B-scan image before any preprocessing. Scale bar: 1 mm for $z$ axis and 10 mm for $x$ axis. B-scan image (b) after applying the thresholding step, (c) showing the result of the connected component analysis, and (d) showing the final segmentation result, with the red line indicating the polynomial fitting curve of the corneal surface. The yellow line indicates the initial segmentation result based on the connected component analysis for comparison purposes.

The subsequent connected component analysis selects the regions with the largest connected elements of the maximal intensity in the image and erases the remaining bright and isolated pixels, retaining the basic shape of the sample boundaries, as illustrated in Fig. 3(c). To further improve this initial segmentation result, which is obtained from low-pass (Gaussian) filtered images, we also high-pass filtered the original B-scan images to obtain tissue boundaries. This result is used to further refine the segmented boundary result. Specifically, we first obtain high-pass filtered images ($F$) by applying a combination of Sobel filter, $S={[\mathrm{1,2},1;\mathrm{0,0},0;-1,-2,-1]}^T$, and Prewitt filter, $P={[\mathrm{1,1},1;\mathrm{0,0},0;-1,-1,-1]}^T$ as follows:

Then, we search the high-pass filtered image, in the vicinity ($\pm 15\text{pixels}$) of the initially segmented surface boundary [yellow line, Fig. 3(d)], for the pixel with the highest gradient (both horizontally and vertically) and use the result to replace the initial segmentation result. The resultant, corrected sample surface segmentation is fitted with a higher-order polynomial curve [red line, Fig. 3(d)]. The polynomial fitting algorithm uses the pixel location to form a Vandermonde matrix ($V$), resulting in a linear system $VP=Z$, which is also

Note that we define the points on the boundary with $(x_i,z_i)(i=1,2,\dots ,m)$, where $m$ is the number of points and $n$ is the order of the polynomial. The coefficients of the polynomial are solved by $P=(V^T{V}^{-1})V^{T}Z$, the matrix multiplication between the Moore–Penrose inverse of $V$ and $Z$.

2.2.2.

Lateral motion-compensation algorithm

The lateral motion-compensation algorithm is based on the preprocessed (segmented) result. We define the first B-scan cornea surface ($z_1$) as the target surface, then move every other B-scan surface ($z_t$, $t$ is the index of B-scan) horizontally to align it with the target surface. For each B-scan, the $n$’th order polynomial fitting curve of the corneal surface can be expressed as

The lateral movement of $t$’th B-scan could be figured out by minimizing cost function

Figure 4 illustrates the lateral motion-compensation mechanism. Figure 4(a) depicts three examples B-scan images (1st, 50th, and 100th) from a volumetric OCT data set of a static ex vivo bovine corneal sample. Figure 4(b) shows the same bovine corneal B-mode images with indicated lateral movement. Figure 4(c) shows the lateral motion-compensated result of (b). The solid blue line in every B-mode image indicates the location of the corresponding corneal surface peak. The dashed blue, red, and green lines in Figs. 4(a)–4(c) are used to illustrate the relative lateral shifts of different B-scans [e.g., $\mathrm{\Delta }{l}_{50}$ and $\mathrm{\Delta }{l}_{100}$ in Fig. 4(b)]. About 50 pixels of both lateral image borders are cropped to remove the empty A-scans generated by the compensation.

Fig. 4

Pictorial explanation of lateral motion compensation by means of an ex vivo bovine cornea. (a) Three static B-scan images serving as a reference. (b) and (c) B-scan images pre- and post-lateral motion compensation, respectively.

2.2.3.

Axial motion-compensation algorithm

The axial motion-compensation algorithm uses the lateral motion-compensated result and is illustrated in Fig. 5. We first compare the B-scan surface detection result with the three reference images and obtain the differences in the axial position at the three intersecting planes of each B-scan. Then, the axial shift for every A-scan within the B-scan is estimated by fitting it with a polynomial and the shift is compensated accordingly in the final images. The motion-free volumetric image can be generated from sequentially-processed B-scans.

Fig. 5

Axial motion-compensation procedure illustrated by means of an ex vivo bovine cornea. (a) Preprocessed and laterally-compensated B-scan image before axial motion compensation, showing the higher-order polynomial fit (blue line), $z(x)$, obtained from preprocessing and the three reference points (yellow dots), $r_i(x_i,y_i,z_i)(i=\mathrm{1,2},3)$, collected from the intersection planes. (b) B-scan image after axial motion compensation, with corresponding final surface fit (red line), $f(x)$.

The three reference points [yellow dots in Fig. 5(a)] are defined as $r_i(x_i,y_i,z_i)(i=1,2,3)$, which we write as $r_i(x_i^r,z_i^r)(i=1,2,3)$, assuming all points in the same B-scan images share the same $y$-axis coordinate. Then, we can deduce the axial motion at the three reference points, $a(x_1^r)$, $a(x_2^r)$, and $a(x_3^r)$ as

The axial motion at every (other) points of $z(x)$, defined as $m(x)$ is estimated by generating a second-order polynomial function based on the three reference points. Finally, we axially shift the location of every A-line in the B-scan to obtain the motion-compensated result [red line in Fig. 5(b)] $f(x)$ as

3.1.1.

Determining order of polynomial for corneal surface fit

With the proposed OCT system being geared toward ophthalmic applications in patients (e.g., to guide ocular surgeries and/or to extract quantitative information from volumetric measurements), our algorithm must work not only on normal corneas but also on pathological ones with abnormal surface topography. We hence used representative clinical OCT images (Avanti, Optovue), including from patients with bullous keratopathy and keratoconus to determine the optimal order of polynomial for the final surface fit during preprocessing of the data. Images were collected retrospectively; however, all subjects provided informed consent to have their images used in research, which was approved by the Institutional Review Board (Patient Protection Committee, Ile-de-France V) and adhered to the tenets of the Declaration of Helsinki. We applied different polynomial fits ranging from second- to sixth-order to the data (image) sets, during step (d) of the preprocessing algorithm [Fig. 2(d)]. Figure 6 shows the results for three typical examples from each group, i.e., for a healthy cornea [Fig. 6(a)], a cornea with bullous keratopathy [Fig. 6(b)], and a cornea with keratoconus [Fig. 6(c)].

Fig. 6

Determining the optimal order of polynomial for the segmented corneal surface fit during the preprocessing step in Fig. 2(d), using clinical SD-OCT (Avanti, Optovue) images of corneas with normal and abnormal surface topology. (a)–(c) The curve fitting with the different order of polynomials (from second to fourth order, as indicated by different colors); (d)–(f) corresponding fitting error, i.e., the RMSE between polynomial fit and manually segmented surface, of the (d) healthy cornea, (e) cornea with bullous keratopathy, and (f) cornea with keratoconus.

We used the fitting error, more precisely the root mean square error (RMSE) representing the difference between the surface fit and the ground truth, with the latter being estimated by manually segmenting the corneal surface, to quantify the results. The resulting fitting error for the three representative corneal images as a function of polynomial order is depicted in Figs. 6(d)–6(f). As expected, the fitting error decreases with increasing polynomial order for the healthy cornea with normal surface topography [Fig. 6(d)]. However, this is not the case for the diseased corneas with abnormal surface pathology [see Figs. 6(e) and 6(f)], where one can note the error starting to increase again for the fifth- and sixth-order polynomials for the cornea with bullous keratopathy [Fig. 6(e)] and keratoconus [Fig. 6(f)] respectively. This may be due to the higher-order polynomials being very sensitive to the small, irregular surface variations of these diseased corneas. The fourth-order polynomial was hence chosen for the final surface fit in step (d) of the preprocessing procedure [see Fig. 2(d)].

3.1.2.

Determining the number of reference planes

As detailed in Sec. 2.2.3, we deduce axial motion associated with every A-line within B-scan images with respect to three reference image planes. To arrive at three as the optimal number of reference planes for axial motion compensation, we tested our algorithm during C-scan imaging (consisting of 512 B-scans) of an ex vivo bovine cornea while inducing typical axial motion of around 1 Hz (such as induced from heartbeat) and changing the number of reference planes from 2 to 4. The reference planes are distributed based on the following rule, which minimizes the compensation error (as discussed in Sec. 3.1.3): suppose there are a total of $t$ planes and $L$ is the length of B-scan images, the $i$’th plane will be located at $(2i-1)L/2t$. The time cost of the polynomial fitting process of the references [described by Eq. (3)] is measured and compared (using the built-in function in C++): $O(n^3)$, in which $n$ is the dimension of the square matrix $V$ in Eq. (3), and it is also the number of references in our case. In our test, first-, second-, and third-order polynomials were generated to estimate axial motion based on the 2, 3, and 4 reference points (planes), respectively. The RMSE between the motion-compensated and static (reference) C-scan image (calculated as the mean value of the absolute difference between the respective A-scan corneal surface positions) was used to quantify the results.

As summarized in Table 1, three references decrease the RMSE by around $4.21\mu \mathrm{m}$ with only an additional 2-ms time cost, compared with two references. Although four reference planes would even further reduce the RMSE by around $0.76\mu \mathrm{m}$, it requires an additional 12 ms for processing and $0.76\mu \mathrm{m}$ is $<1\text{pixel}$ and the resolution of our imaging system. Hence, the three reference-plane approaches provide an optimum compromise between performance and time cost for our overall scanning speed and considered frequency range.

Table 1

RMSE and time cost for C-scan imaging with different numbers of B-scan references.

3.1.3.

Determining the reference plane positions, $\alpha$

Theory. The compensation result also depends on the position of the reference images and the distance between the adjacent reference images, $d$, and thus our selection of $\alpha$, which is defined as $d/L$, where $L$ is the width of the raw B-scan image [see Fig. 1(b)]. Since we choose the middle reference to be fixed at the center of the $x$-direction, $\alpha$ can range from 0 (i.e., one reference plane at the center) to 0.5 (i.e., two references at the ends of the image).

This is schematically shown in Figs. 7(a)–7(c), displaying the location of the three references ($R1$, $R2$, and $R3$; dotted lines) in a B-scan image, with $\alpha$ set as 0, 1/3, and 0.5, respectively.

Fig. 7

Analysis for the selection of. (a)–(c) Schematically showing the location of the three references (black dotted lines) in a B-scan image, with $\alpha$ set at 0, 1/3, and 0.5, respectively. (d)–(f) The corresponding accuracy distribution of A-scans within the B-scan. Black dotted lines indicate the location of the three references ($R1$, $R2$, and $R3$).

The quantitative analysis model is based on the hypothesis that the motion compensation accuracy based on references is a Gaussian-like function of the number of A-scans within the B-scan. The fitting error within every B-scan image is considered random, so the accuracy of references used to compensate for the motion could be estimated as a Gaussian function (i.e., a widely used noise model to describe estimation-related error²⁶), and our experimental result agrees with this hypothesis.

The accuracy function may be written as

We assumed $m=1023$ to be the maximum number of A-scans. All three reference functions should have the same $\sigma$ because all of them have the same compensation capability. Therefore, we have $\sigma ={\sigma }_i(i=1,2,3)$.

We define another accuracy function $Y(x)$ to determine the total accuracy of the motion compensation based on more than one reference. $Y(x)$ is a combination of several $R_i(x)$ and should follow two rules:

These two rules ensure that $Y(x)$ also ranges from [0,1]. Thus, we construct $Y(x)$ as

The number 1 in the denominator of distribution function $(1/(|x-{\mu }_i|)+1)$ is used to avoid zero. Figures 7(d)–7(f) show the accuracy functions $R_1(x)$, $R_2(\mathrm{x})$, $R_3(x)$, and $Y(x)$ with $\alpha$ set to 0, 1/3 and 0.5, respectively. Among the three examples, we obtain the best accuracy result with $\alpha =1/3$. The corresponding minimum accuracy of 0.33 represents the maximum error when the speed of the motion is approximately the same as the scanning rate of the reference images. However, since the reference images are obtained around 100 Hz, which is much faster than the anticipated sample motion, ranging up to 5 Hz, we expect our overall result to be highly accurate.

We also define an error function based on the accuracy function above to enable comparison with experimental data

Combining Eqs. (8) and (10) and assuming $\sigma$ as a constant, we can write $E({\mu }_i,\sigma )$ as a function of $\alpha$

Experiment. To optimize $\alpha$ experimentally, we assessed the performance of the motion-compensation algorithm by performing several imaging trials while varying the distance between adjacent reference planes. The RMSE representing the difference between respective A-scan corneal surface positions of the motion-compensated and static (reference) C-scan image was used to quantify the results. Figure 8(a) shows the RMSE of the compensation result for around 1-Hz motion using two bovine corneal samples as a function of $\alpha$. Two trials were performed for each sample. Figure 8(b) shows the RMSE result for different motion frequencies (0.4 to 5.2 Hz) using the corneal phantom as a function of $\alpha$. Both experiments indicate that $\alpha$ around 0.33 achieves the best results (i.e., minimal RMSE) for all samples and frequencies used in the study. Note that, while the slower motion frequencies yielded smaller compensation errors (i.e., an RMSE of 2.34 for 1 Hz, and an RMSE of 1.26 for 0.4 Hz), even at higher frequencies (2.4 through 5.2 Hz), with $\alpha$ set to 0.33, our algorithm was capable of compensation with an error of $<4.2\mu \mathrm{m}$. Note that although the axial resolution of our system is $6\mu \mathrm{m}$ in air and $4.5\mu \mathrm{m}$ in corneal tissue, our ability to detect the peak of the point spread function is limited by the pixel size, which is $3.6\mu \mathrm{m}/\text{pixel}$ (for corneal tissue). Hence, $4.2\mu \mathrm{m}$ corresponds to an overall compensation error of about 1 pixel.

Fig. 8

RMSE between motion-compensated and static reference images as a function of $\alpha$ for (a) two different ex vivo bovine cornea samples with two different trials with 1 Hz axial motion and (b) a phantom result with motion frequencies from 0.4 Hz to 5.2 Hz.

Comparison. To further test the performance of our model, we calculated the error function with different $\sigma$ [using Eq. (10)] and compared it with the experimental data, where the difference is the fitting error defined as

As illustrated in Fig. 9(a), the fitting error ($W$) achieves a minimum at $\sigma =123$. Figure 9(b) depicts the result, at minimal $\sigma =123$, as a function of $\alpha$ (step size 10/1024), for the experimentally obtained motion compensation result and the numerical result for a motion frequency of 5.2 Hz; note that the data has been normalized. The average absolute error difference between the theory and the experiment is 0.11. Both theory and experiment achieve a minimum error around $\alpha =0.338$. Hence, we conclude that the mechanism of our motion compensation method is successfully described by our theoretical model.

Fig. 9

(a) Fitting error ($W$) between our theory and experimental data for a motion at 5.2 Hz as a function of $\sigma$. (b) Comparison between the experimental motion compensation result (black) and numerical result using the theoretical model (red) as a function of $\alpha$, for minimal $\sigma =123$.

The result above shows that the evenly-distributed reference planes can minimize the compensation error, which can be generalized to the distribution rule of reference planes we mentioned in Sec. 3.1.2.

3.2.1.

Fourth-order polynomial for corneal surface fit

To validate the selection of the fourth-order as the optimal polynomial for the corneal surface fit, we assessed the performance of the algorithm experimentally. Specifically, two imaging trials using two ex vivo bovine corneal samples while inducing quasisinusoidal axial motion around 1 Hz were performed. The resulting images were fitted with the different order polynomials and compared. As previously done, the RMSE was used to quantify the results. The RMSE was deduced from measuring the difference between the reference and compensated C-scan image of the sample. The results are illustrated in Fig. 10 and show that fourth-order fitting achieved the best accuracy, i.e., the minimum RMSE. Considering this and our preliminary result using clinical images of human corneas with normal and abnormal surface topology, the fourth-order polynomial was confirmed as optimal surface fit in step (d) of the preprocessing procedure [Fig. 2(d)].

Fig. 10

RMSE for four compensation studies using two cornea samples while inducing axial motion of 1Hz as a function of polynomial orders from 2 to 6.

3.2.2.

Proof-of-concept demonstration of motion compensation

Using the optimized parameters, we performed a series of volumetric imaging trials using corneal phantom (thin plastic shell) and tissue samples (ex vivo bovine cornea), induced with artificially generated motions, as a proof-of-concept demonstration of the proposed motion-compensation method using our in house-built 3D-OCT system.

Axial motion. Figures 11(b), 11(e), 11(h), and 11(k) depict the volumetric (C-scan) images of the corneal phantom during induced axial motion [by means of the VTS, Fig. 1(a)] around 0.4, 2.4, 3.6, and 5.2 Hz, respectively, clearly showing the motion artifacts. Figures 11(a), 11(d), 11(g), and 11(j) are the corresponding static C-scan images, without any induced motion, that were used as references, and Figs. 11(c), 11(f), 11(i), and 11(l), depict the respective motion-compensated C-scan images using the proposed method. As expected, our algorithm successfully corrects the motion artifacts in these four motion frequency settings, with the smooth spherical surface of the phantom being visibly recovered throughout the entire volume of the sample. The RMSE, as quantitative metric of the results, is 1.3, 3.5, 3.7, and 4.2 μm for the axial motions of 0.4, 2.4, 3.6, and 5.2 Hz, respectively [see Fig. 8(b)].

Fig. 11

Axial motion-compensation result during volumetric (C-scan) imaging of a corneal phantom (thin plastic shell). Panels (a), (d), (g), and (j) are the static reference images of the phantom. The reference planes are indicated in red in (a). Panels (b), (e), (h), and (k) correspond to the sample motion around 0.4, 2.1, 3.6, and 5.2 Hz, respectively. Panels (c), (f) , (i), and (l) are reconstructed after using the proposed motion compensation method.

An example of the axial motion-compensation result during volumetric imaging of an ex vivo bovine cornea is illustrated in Fig. 12. The average amplitude and frequency of the induced corneal motion were around $700\mu \mathrm{m}$ and 1 Hz, respectively. Figure 12(a) depicts the C-scan reference image of the cornea without any induced axial motion, whereas Figs. 12(b) and 12(c) are the C-scan images of the “moving” cornea before and after motion compensation. As expected, the motion-induced distortion of the corneal surface, clearly visible in Fig. 12(b), is successfully compensated using the proposed method [see Fig. 12(c)]. This example corresponds to trial 1 of sample 1 in Fig. 8(a), with an RMSE of $1.35\mu \mathrm{m}$ (for $\alpha =0.33$). For direct visual comparison, the slow-axis cross-section showing the difference between reference and compensated result is also shown in the inset of Fig. 12(c).

Fig. 12

Axial motion-compensation result during volumetric (C-scan) imaging of an ex vivo bovine cornea. Panel (a) is the static reference of the cornea; panel (b) is the cornea with simulated axial motion of around 1 Hz. Panel (c) is the reconstructed corneal C-scan image after motion compensation. The red line in (a) and cyan line in (c) represent the corneal cross-section location of reference and compensation result, respectively. Both are also shown in the inset of (c), which provides a direct comparison between reference and postcompensation results.

3D motion. Figure 13 shows an example of the 3D motion-compensation result during C-scan (volume) “free-hand” imaging of a phantom [plastic shell; Figs. 13(a)–13(c)] and ex vivo bovine corneal sample [Figs. 13(d)–13(f)], i.e., with the samples positioned in the palm of the hand, to more realistically simulate the involuntary motion of living tissue occurring in both lateral and axial directions. Figure 13(a) and 13(c) depict the static (reference) C-scan image of the phantom and bovine cornea, respectively, without any motion; Figs. 13(b) and 13(e), and Figs. 13(c) and 13(f) are the C-scan images of the phantom and held corneas, before and after motion compensation, respectively. The marked central line on the corneal surface is a manually produced scar (by scratching the surface) to highlight the occurrence of lateral and axial motion during free-hand imaging [see Figs. 3(b) and 3(e)]. Note that the marked line appears straighter again after compensation [see Figs. 3(c) and 3(f)]. The RMSE of the phantom and bovine cornea result is 4.61 and $3.96\mu \mathrm{m}$, respectively.

Fig. 13

3D motion-compensation result of a phantom (top) and ex vivo bovine (bottom) cornea during C-scan (volumetric) free-hand imaging, i.e., with the samples positioned in the palm of a hand to simulate involuntary tissue motion in both lateral and axial directions. Panels (a) and (d) are the static reference of the cornea without any simulated motion; panels (b) and (e) are the cornea with free-hand motion. Panels (c) and (f) are the reconstructed corneal image after 3D-motion compensation.

We also analyzed the lateral and axial motion that was induced during free-hand imaging of the corneal samples. The amplitudes of axial and lateral motion occurring during free-hand imaging of the ex vivo bovine cornea from Fig. 13 are depicted in Figs. 14(a) and 14(b), respectively; the corresponding frequencies are shown in Fig. 14(c). The mean value of axial and lateral motion between adjacent B-scan images is around 14.18 and $165\mu \mathrm{m}$, respectively. Also note that most energy of motion, in either direction, occurs in the range from 0 to 5 Hz and that motion above 15 Hz is negligible.

Fig. 14

Axial and lateral motion analysis during free-hand volumetric imaging of an ex vivo bovine cornea. Panels (a) and (b) are the amplitude of axial and lateral tissue motion, respectively; panel (c) is the corresponding frequency information.