1. Introduction

Ophthalmologists rely on accurate corneal topography for surgical planning and diagnostics, e.g. to detect pathological deformations, such as keratoconus, at an early stage. Corneal topography provides two-dimensional maps of the cornea’s curvature and refractive power. Placido disc reflection, scanning slit and Scheimpflug photography are established modalities for corneal topography. While Optical Coherence Tomography (OCT) is commonly used for cross-sectional or volumetric imaging with high resolution, its use for corneal topography is still limited. In contrast to conventional modalities, OCT potentially enables a simultaneous measurement of all optically relevant structures of the eye, including, but not limited to, corneal topography. Whereas traditional photography-based principles enable acquisition of 2-D or 3-D geometrical information of the cornea in an instant, beam-scanning OCT relies on the sequential gathering of one-dimensional axial depth profiles (A-scans) to obtain two-dimensional scans (B-scans) or volumetric scans. For OCT-based corneal topography, the cornea is scanned at discrete positions around the corneal vertex. This sequential scanning takes time and makes the measurement susceptible to eye motion. The acquisition time could be minimised by using high-speed OCT systems, yet these systems are costly and the high acquisition rate sacrifices signal quality. The need for lateral scanning could be entirely removed by using parallel OCT instead of beam-scanning OCT, but these systems have other disadvantages, such as the need for mechanical depth scanning or more complex optical setups [1,2]. Real-time tracking of the corneal position during acquisition is another rather hardware-demanding and expensive approach. Therefore, common approaches rely on tailored scanning schemes and algorithms to enable motion compensation after the acquisition.

The feasibility of motion compensation is strongly linked to the used scan pattern. Recently, Wagner et al. [3] proposed golden angle based spiral scanning for robust corneal tomography. This scanning scheme enables the elimination of abrupt and large movement from the measurement without sacrificing the reliability of the reconstruction. However, the problem of continuous motion during scanning was not addressed. Scan patterns more commonly used for corneal topography consist of several meridional B-scans centred around a common central point which is usually aligned to the apex of the cornea [4–7]. These so-called radial scan patterns simplify axial motion compensation because all meridional scans share the apex as a common point, allowing an alignment of the individual B-scans. McNabb et al. [8,9] use a special radial scan pattern for their Distributed Scanning OCT method, consisting of 20 meridional B-scans with 500 A-scans each. Spatially adjacent A-scans are not gathered sequentially in time but distributed over five sub-sampled passages. Thus, the scan pattern consists of five connected scans with 20 radial B-scans containing 100 A-scans. Each B-scan takes 10 ms and is assumed to be free of motion. Motion-compensated meridional B-scans are then reconstructed from the five sub-sampled B-scans along the given meridian. Raster scanning is another, more traditional scanning type. Although raster scanning is not commonly used for corneal topography, motion compensation methods are applied to generate distortion-free volumes of the retina or anterior segment. In classical raster scanning, several parallel B-scans are gathered, covering a rectangular area. This results in a fast and a slow scan axis, one perpendicular to the other. While B-scans along the fast scan axis are considered motion free [10], B-scans along the slow axis are prone to motion. Although visually appealing volumes can be created by aligning the fast-axis B-scans to each other, geometrical information along the slow scan axis will be lost. To overcome this problem, Kraus et al. [11] gathered additional B-scans along the slow scan axis to align the fast-axis B-scans to these additional B-scans. Recently, Chen et al. [12] presented a Lissajous scan pattern for three-dimensional imaging of the retina. Because the scan pattern frequently overlaps with itself, they were able to compensate motion in three dimensions.

We present a novel pipeline for motion-compensated corneal topography by OCT, including the following steps (see Fig. 1): (1) Scanning of the cornea with a tailored scan pattern, (2) chunk-wise segmentation of the anterior corneal interface, (3) continuous motion compensation in three dimensions, (4) outlier removal, (5) regularised Zernike reconstruction and (6) generation of the axial curvature map. Our scan pattern provides fast coverage and a high number of intersections to aid motion compensation. The motion compensation features continuous, model-based motion detection in three dimensions with high temporal resolution. In contrast to common methods [8,11,12], no part of the scan is assumed to be motion-free. Instead, our method features continuous motion compensation in three dimensions. We model the motion by a piecewise polynomial spline and formulate an optimisation problem to determine the motion parameters together with a parametric description of the anterior corneal surface. The regularisation of the reconstruction addresses instabilities in the curvature space of standard Zernike reconstruction.

 

Fig. 1. The different steps in our pipeline.

Download Full Size | PDF

2. Methods

2.1 Scan system

We use a prototype swept-source OCT system optimised for simultaneous OCT imaging of the anterior eye segment and measurement of axial eye length in a single scan range. Figure 2 shows a schematic of the system with its main components. The system employs a telecentric scan geometry. To achieve sufficient collection efficiency over the entire eye length, we use a small numerical aperture of approximately NA $\approx 0.015$. The focal plane is 5 mm behind the corneal vertex (neglecting the refractive effect of the eye). As a result, the lateral resolution at the anterior cornea is approximately 80 µm (full width at half maximum).

 

Fig. 2. Schematic of the swept-source OCT system with its components. LS: Laser Source; PC1, PC2: Polarisation controllers; FC: Fibre coupler; CL: Collimating lens; MS: 2-D MEMS beam scanner; MO: Measurement objective; D: Detector; FM: Fibre mirror.

Download Full Size | PDF

The laser source (LC, Axsun Technologies) sweep is centred at 1060 nm with a sweep range of 30 nm and a sweep rate of 30 kHz. The coherence length of the laser is 75 mm.

We comply with the requirements specified in ISO 15004-2:2007 (Group 2) [13] to limit the optical power. With a sample arm power of 1.9 mW, the OCT system has a sensitivity of 104.4 dB at the focal plane, 10 mm from zero delay. The total axial measurement range is 50 mm in air, and the axial resolution is 40 µm in air (full width at half maximum). We employ a 2D MEMS scanner (MS, Hamamatsu Photonics) for lateral scanning. Scan patterns with reduced bandwidth are used to avoid excitation of MEMS resonances around 500 Hz and 1 kHz. Since the scanner is operated in open-loop (without position encoders), careful calibration of the scanner control signal and characterisation of the resulting scan positions are required, as described in the following sections.

2.2 Scan pattern

As we mentioned in the introduction, it is common to use intersecting scan patterns to enable motion compensation: Additional B-scans along the slow axis of raster scans are used to align the B-scans to match at the intersections. Radial scan patterns facilitate motion compensation because the individual B-scans intersect at the vertex. The distributed scan pattern of McNabb et al. increases temporal and spatial linkage by gathering spatially adjacent scans distributed over five sub-sampled passages – also to enable motion compensation. The same principle applies to the Lissajous scan pattern used by Chen et al. [12].

We developed a new scan pattern that combines fast coverage with a high number of intersections while keeping the demands on the dynamics of the scanner low. The high number of intersections results in high temporal and spatial linkage in the scan. Because of this linkage, motion leads to apparent inconsistencies in the measurement, which enables assessment of the motion. The scan pattern trajectory is defined as

 

Fig. 3. (a) One frame of our scan pattern and (b) its projections onto the $x$ and $y$ axis. The color coding of the scan points indicate the point in time and is consistent between both subfigures.

Download Full Size | PDF

We calibrated and characterised the scan pattern based on scans of a flat mirror as described in Section 2.3. The mirror is placed at a scan depth of 5 mm, which corresponds to the depth of the corneal vertex when scanning the cornea. The mean absolute deviation between the calibrated scan points and the specified scan points was 37 µm with a maximum of 73 µm. The deviation between specified scan points and calibrated scan points are mainly due to distortions caused by scan optics. However, the calibration method fully characterises the distorted pattern. The uncertainty of the calibrated scan points due to scan-to-scan variations and uncertainties of the calibration method is around 5 µm.

2.3 Calibration and characterisation of the scan pattern

For the applicability of the presented motion correction and reconstruction methods, the calibration method is not crucial, and other approaches can be used instead [5,14]. However, since the accuracy of lateral beam positioning can limit the accuracy of corneal curvature measurement by OCT, the calibration and characterisation of the scan pattern require special attention, in particular in an open-loop scan system as used in this work. Since to our knowledge, the method applied here has not been discussed in literature before, we provide a brief outline. The method is based on a measurement of scan point coordinates by the OCT system itself (“interferometric scan position measurement”) and consists of three steps: (1) Characterisation of the two-dimensional MEMS transfer function, (2) calculation of the scanner control signal for the scan pattern, (3) characterisation of the resulting scan positions.

2.3.1 Interferometric scan position measurement

We place a flat mirror on a gimbal mount in front of the system, centered with respect to the scan range of the beam scanner. To determine the lateral positions of a scan point, we measure the distance $z_{1,2,3,4}$ to the surface at that scan point by OCT, at four different tilt angles of the surface normal with respect to the beam axis, $(\Theta _x, \Theta _y)_{1,2,3,4} = (+\alpha , 0), (-\alpha , 0), (0, +\alpha ), (0, -\alpha )$. The $z$-axis is parallel to the beam axis, and $x$ and $y$ are orthogonal to the beam axis and with respect to each other. The lateral position of the scan point is then given by $x = ( z_2 - z_1) / (2 \tan \alpha )$, and $y = ( z_4 - z_3) / (2 \tan \alpha )$. The axial offset of the scan point is given by $z_{\mathrm {ofs}} = ( z_2 + z_1) / 2$, or $z_{\mathrm {ofs}} = ( z_4 + z_3) / 2$. To take the frequency dependence of the scanner deflection into account, we measure a complete period of the scan pattern (points $i=1\ldots N$) for each mirror orientation and calculate the resulting trajectories for $x_{1\ldots N}$, $y_{1\ldots N}$, and $z_{\mathrm {ofs,1\ldots N}}$, based on the $4 N$ distances $(z_{1,2,3,4})_{1\ldots N}$. The necessary precision of the distance measurement is achieved by Gaussian peak fitting to each A-scan.

2.3.2 Characterisation of the two-dimensional MEMS transfer function

For each scanner axis individually, we measure the beam trajectory resulting from a small amplitude sinusoidal input function using the method described in Section 2.3.1, yielding the response amplitudes in $x$ and $y$-axis. This way, we can build a $2\times 2$ matrix mapping a control signal ($c_1$, $c_2$) to a beam position ($x$, $y$) in the frequency domain. We repeat this for a set of frequencies between 25 and 300 Hz, and fit each matrix element by a third-order polynomial of the frequency, yielding the two-dimensional MEMS transfer function.

2.3.3 Calculation of the scanner control signal

We convert the scan pattern into the frequency domain using a Fourier transformation (the scan pattern is periodic). We then apply the inverted two-dimensional MEMS transfer function and convert the resulting control signal back into the time domain using an inverse Fourier transformation.

2.3.4 Characterisation of the resulting scan positions

We apply the control signal obtained according to Section 2.3.3 to the scanner, and characterise the resulting scan point coordinates according to Section 2.3.1. For the OCT system discussed in this work, the standard deviation of repeated characterisations of the scan pattern per point is typically 5 µm (averaged over all scan points). The mean absolute deviation between the characterised and specified scan point coordinates was 37 µm with a maximum of 73 µm. This deviation is mainly due to distortions caused by scan optics. It can be properly taken into account when evaluating corneal curvatures by using the characterised coordinates instead of the originally specified coordinates so that it does not harm the accuracy of the corneal topography measurement.

2.4 Segmentation

We use a graph-based segmentation, similar to the one described by Wagner et al. [15], to segment the anterior corneal surface in the scans. The key feature is a model-driven, three-dimensional regularisation of the region of interest (ROI), slope and curvature. For the segmentation, we divide each frame of 16384 A-scans into two chunks. Hence, we split the complete scan of four frames into eight chunks of 8192 A-scans each. For each chunk, the A-scans containing eyelid are identified and excluded from the segmentation (cf. Fig. 4). The automatic eyelid classification uses the centre of mass of the A-scan intensities (first moment of the A-scan vector divided by the sum of the A-scan vector) as the main classification feature. The anterior corneal surface is then segmented chunk-wise by using a regularisation model that is continuously updated between the chunks.

 

Fig. 4. Section of the oversampled gradient image with masked out eye lid (reddish overlay), the ROI (blue) and the fine segmentation of the anterior corneal surface (green).

Download Full Size | PDF

In general, we use the segmentation of each chunk to update the regularisation model for the segmentation of the next chunk. For the first chunk, we initially segment the corneal stroma using an initial model of the cornea for regularisation. Based on the stroma segmentation, the model for the segmentation of the anterior surface is generated. Wagner et al. [15] described the initial model, the preprocessing and the segmentation procedure in more detail. Each segmentation consists of a coarse and a fine segmentation. The coarse segmentation is used to create a refined model for the fine segmentation. It relies on the downsampled (by a factor of 4) gradient image, whereas the fine segmentation uses an oversampled (by a factor of 2) gradient image, providing a segmentation value for every 8th A-scan.

After the segmentation of the anterior corneal surface in a chunk, the segmentation is checked for plausibility. Implausible segmentation can be caused by blinking of the patient, saccades during the chunk or extensive motion. For the plausibility check, a 5th-order Zernike reconstruction is performed on the segmented points of the chunk, and the mean coefficients standard error [15] is used for classification. The threshold for classification was set to 2.5×10^-6 m. Please note that this threshold is only valid for the given chunk size and the expected number of segmented points. If classified as implausible, the segmentation of this chunk is dumped, and the regularisation model is not updated. In this case, the old regularisation model is used for the segmentation of the next chunk.

Although we did not explicitly implement detection of big saccades or extensive motion, the combination of the model-based ROI and exclusion of implausible chunks results in removal of chunks affected by big saccaded or extensive motion: If big saccades or extensive motion are present in a chunk, the cornea will be outside the ROI. This leads to segmentation in the noise, which further results in a high mean coefficients standard error and exclusion of the chunk. Chunks containing distortions not resulting in failed segmentation but with a negative effect on the reconstruction are removed during outlier removal (see Sec. 2.6).

Because the signal of the phantom is different from the signal of the cornea, we use a modified segmentation scheme for phantoms: First, no initial stroma segmentation is performed. Instead, the phantom surface is segmented directly using the initial model. Second, the segmentation does not rely on a gradient image. Instead, we rely on the same preprocessing that was used for phantoms by Wagner et al. [15].

2.5 Motion compensation

Our motion compensation uses the segmented points to estimate and compensate motion-induced displacement in all three axes (see Fig. 5) simultaneously, while we assume the rotations to be negligible. We employ a continuous optimisation scheme to reproduce the position of the segmented points by Zernike polynomials and three-dimensional piecewise polynomial displacement. The optimisation, therefore, yields a polynomial description of the displacement in each chunk. The displacement is forced to be continuous and continuous in its first derivative between adjacent chunks by applying the corresponding regularisation.

 

Fig. 5. The movement of the eye and its axes.

Download Full Size | PDF

2.5.1 Optimisation problem

The point cloud from the segmented anterior corneal surface is the input of the optimisation. For each segmented point, we have the coordinate $P_n = \begin {pmatrix} x_n & y_n & z_n \end {pmatrix}^T$ and the time $t_n$ for $n = 0, 1, \ldots , N^t-1$, where $N^t$ is the total number of points. As for the segmentation, we split the measurement into $Z=8$ intervals of equal duration, which results in 8 point chunks. Because the segmentation can have gaps due to removed eyelids, the number of points in the chunks varies. Empty point chunks, caused by implausible segmentation (cf. 2.4), are excluded from the following steps, resulting in $Z^* \leq Z$ valid point chunks. The piecewise polynomial displacement is defined by

We determine the coefficients of the displacement polynomials by solving the linear least-squares problem

2.6 Outlier removal

In addition to the removal of implausible chunks during segmentation, we remove outlying chunks in a step between motion compensation and reconstruction. The outlier-removal targets chunks that suffer from saccades or other distortions that were not compensated well enough by the motion compensation and are too small to be detected by the plausibility check during segmentation. Because microsaccades occur at an average rate of 1 to 2 per second, we expect microsaccades in 2 to 4 chunks of the 8 chunks.

The algorithm iteratively removes chunks that increase the coefficients standard error[15] of the overall Zernike reconstruction. Given a base set of $N$ segmented and motion-compensated chunks, $N$ different subsets containing $N-1$ chunks are generated by excluding single chunks. The coefficients standard error of the 7^th-order Zernike reconstruction is determined for each subset. If the smallest coefficients standard error is smaller than the coefficients standard error of the overall reconstruction, the procedure is repeated with the corresponding subset as the new base set. If not, no outlier is present, and the base set is passed to the reconstruction.

2.7 Reconstruction

We use Zernike polynomials to reconstruct the anterior corneal surface from the motion-compensated points and to determine the radial curvature of the surface at points relative to the topographers axis – which is typically aligned to the visual axis. The Zernike polynomials are defined by

We sampled $Z$ at the individual scan positions $x_i, y_i$ and plot the resulting reconstructed axial distances $Z_i=Z(x_i,y_i)$ in the scan (see Fig. 6, yellow line). As the A-scans in the scan are displayed as recorded, without any motion correction, the reconstructed axial distances do not coincide with the scanned anterior corneal surface. To illustrate the accuracy of our reconstruction, the compensated motion $f$ has to be incorporated by $Z_i^*=Z\left (x_i+f_x(t_i),y_i+f_y(t_i)\right ) + f_z(i)$ (cyan line). The match of the cyan line with the scanned anterior corneal surface indicates accurate motion compensation and reconstruction. The difference between the yellow and cyan line illustrates the compensated motion.

 

Fig. 6. Section of a corneal scan as acquired (no motion correction applied). The reconstruction of the anterior corneal surface is delineated by sampling the reconstructed corneal model $Z$ at the scan positions $x_i, y_i$. Yellow: Sampling of the static model, representing the axial distances to the virtual motion-free anterior corneal surface ($Z_i$). Cyan: Sampling of the model incorporating the motion profile $f$ determined during motion compensation, recreating the axial distances, as acquired ($Z_i^*$).

Download Full Size | PDF

The axial curvature $K_{\mathrm {a}}$, defined by the ISO norm for corneal topographers [17], is the reciprocal of the distance from a point on the surface to the corneal topographer axis along the corneal meridian normal. It is defined by

2.8 Validation

2.8.1 Phantom validation

We used a toric acrylic glass (PMMA) phantom with a steep (flat) radius of 7.8 mm (8.0 mm) to evaluate the performance of the motion compensation. The phantom was manufactured by Sumipro B.V., Netherlands. The shape of the phantom represents a typical eye with astigmatism. The preprocessing for the segmentation was adapted for the segmentation of the phantom surface (see Sec. 2.4). We applied motion, simulating typical eye displacements to the phantom by using a piezoelectric translation stage. It consists of axial (along the $z$ axis, which corresponds to the topographer axis) and lateral motion trajectories (in the $x$,$y$ plane). The applied axial motion was extracted from a prior in vivo measurement with our system where we scanned the corneal vertex over 1.75 seconds using a static scan pattern. Because we aligned the scanning to the corneal vertex and the corneal slope is minimal around the vertex, we assumed the change of the corneal position in the scan to be caused by axial eye movement only. The motion applied in the $x$,$y$ plane represents drift and microsaccades. The drift was parameterised with amplitudes between 3 µm and 30 µm and velocities between 5 µm and 30 µm per second. Saccades were simulated in intervals between $\frac {1}{3}$ and 1 s. The saccades were parameterised with amplitudes between 9 µm and 150 µm and velocities between 1.4 mm and 14 mm per second.

We acquired 10 measurements from the moving phantom and 5 measurements from the same phantom without any movement.

2.8.2 In-vivo validation

We evaluated the in vivo repeatability and equivalence of our topography by measurements of 48 healthy eyes. The topography for both eyes of 24 subjects was determined by our method and by two reference topographers, a Scheimpflug imaging tomographer (Pentacam HR, Oculus) and a Placido Disk topographer (Atlas 9000, Carl Zeiss Meditec). Each subject was measured twice with each device, where the subject setup was repeated before each measurement. Thus, 48 eyes were measured twice with three devices, resulting in a total of 288 acquired topographies.

We excluded three subjects entirely from the evaluation. For these subjects, the Atlas indicated low measurement quality and our segmentation showed a deviation caused by an unusual corneal signal. We additionally excluded one subject from the Pentacam measurements and three measurements from the Atlas measurements because the particular device indicated a low measurement quality.

Based on the remaining topographies, we evaluated the repeatability of each device and the mutual equivalence between the devices, as explained in the following sections.

The measurements were conducted by the Eye Clinic of the University Hospital Basel with approval of the local ethics committee. Following the tenets of the Declaration of Helsinki, the measurements were performed with the full understanding and written consent of the volunteers.

2.8.3 Calculation of the difference between topographies

For the calculation of the difference between two topographies, we first ensure that the vertices of the topographies are aligned. As the maps provided by the Pentacam and the Atlas are already aligned to the vertex, we only had to align our map. The position of the vertex was determined from the Zernike coefficients as described in [19]. Figure 7 shows two aligned topography maps and the mapped differences. The ISO norm for Corneal topographers [17] defines two measures for the difference between two topographies: the difference mean and the standard deviation of the differences. To compare two topographies $i$ and $j$, both topography maps are sampled at a defined set of positions $k$ and the weighted differences for these positions are calculated by

 

Fig. 7. (b) Axial curvature map of our OCT-based topography. (c) Axial curvature map of the Pentacam. (d) The mapped differences between the two topographies. The axial curvature maps are color-coded with increments of 0.5 D between the shades of color (a). The difference map is color-coded with increments of 0.2 D between the shades of color (e). The pentacam map is cropped to a square area with 8 mm lateral extent.

Download Full Size | PDF

We interpret the standard deviation of the differences as a measure for the difference of the aberrations in two topographies. In total, the calculation of the difference between two topographies yields six numbers: the difference mean and standard deviation of the differences for three zones.

2.8.4 Calculation of the repeatability and equivalence

In vivo repeatability For the assessment of in vivo repeatability, we calculated the difference between the first and second topography of each eye and device. Each difference includes the difference mean and standard deviation of the differences for the three zones.

Based on these differences (per eye), we calculate two measures for the repeatability of each zone (per ensemble): (1) the standard deviation of the difference mean as a measure for the repeatability of the mean, and (2) the mean of the standard deviation of the differences as a measure for the repeatability of the aberrations. We use the mean as a measure for the repeatability of the aberrations because the standard deviation of the differences is always positive.

Phantom repeatability To assess the repeatability of the topography from the 10 measurements of the moving phantom, 10 pairs were generated, and the topography difference was calculated for each pair. Given the 10 topographies $T_i$ with $i=0\ldots 9$, the differences $T_j-T_{j+1}$ for $j=0\ldots 8$ and $T_9-T_{0}$ were determined. The measures for the repeatability of the mean and the repeatability of the aberrations were calculated analogously to the in vivo repeatability. The same procedure was applied to the assessment of the repeatability of the 5 topographies of the static phantom.

In vivo equivalence For the assessment of the equivalence of two devices, we calculated two differences for each eye: the difference between the first topography of each device and the difference between the second topography of each device. Each difference includes the difference mean and standard deviation of the differences for the three zones.

Based on these differences (per eye), we calculate two measures for the equivalence of each zone (per ensemble): (1) The mean of the difference mean as a first measure for the equivalence of the mean, (2) the standard deviation of the difference mean as a second measure for the equivalence of the mean, and (3) the mean of the standard deviation of the differences as a measure for the equivalence of the aberrations.

2.8.5 Adjustment of the topography for equivalence of the mean

We observed systematic differences between the OCT-based delineation of the anterior corneal surface and the Placido disc reflection. These differences were expected in advance due to the different measurement modalities. Thus, we performed a simple adjustment of our topography to the topography provided by the Atlas. The adjustment parameters were determined based on additional in vivo measurements, acquired analogously to the measurements for in vivo validation. The adjustment was applied to $z$-coordinates of the segmented point by a conic correction. For the correction, we first fit a conic section of revolution to the motion-compensated point cloud [18]. The parameters $r$ and $p$ of the conic section of revolution are adjusted by $R^* = R + \Delta R$ and $p^*= p + \Delta p$, resulting in an adjusted conic section of revolution. The difference between the $z$-position of the adjusted and the original conic section of revolution is determined at the lateral position of each point and added to the $z$-coordinate of the point. The reported results are based on an adjustment of $\Delta R = 0.166 \times 10^{-3}$m and $\Delta p = 0.19$,

3. Results

Effect of the regularisation Figure 8 illustrates the effect of the global regularisation weight on the in vivo repeatability of our OCT-based topography and the equivalence with the reference topographers. While the repeatability steadily improves with increasing regularisation, the equivalence with the Altas in the central zone becomes worse after an initial improvement. Based on these characteristics, we have chosen a regularisation weight of $\lambda = 0.004$ for the following evaluation.

 

Fig. 8. Influence of the regularisation on (a) the repeatability of the aberrations and (b) the equivalence of the aberrations with Atlas (solid) and Pentacam (dashed), for the central zone (red), middle zone (green) and outer zone (blue).

Download Full Size | PDF

Phantom repeatability and equivalence Figure 9 illustrates the effect of the motion compensation on the repeatability of the moved phantom topography, in comparison with the repeatability of the static phantom topography. As can be seen, the motion compensation improves the repeatability by orders of magnitude and brings it to the level of the repeatability of static phantom topography. Figure 10 provides an example for the detected and compensated motion in one of the measurements, showing good correlation. We zoomed in a typical response of the motion compensation to a simulated saccade. Although the motion compensation is able to correct the shift induced by the saccade, the speed of the motion compensation response is limited to around 0.1 s, which corresponds to a maximal frequency of the lateral motion compensation of around 10 Hz. Figure 11 illustrates the difference between the theoretical topography of the phantom, the topography of the moved phantom with and without motion compensation, and the deviation between the motion-compensated topography and the theoretical topography. Table 1 lists the corresponding difference measures (see Sec. 2.8.4). As can be seen, the difference mean in the middle and outer zones is low for the topography without motion compensation although the topography map is quite different. This confirms that the measure standard deviation of the differences is needed to quantify the difference between topographies as the other measure, difference mean, is blind for irregular differences (aberrations).

 

Fig. 9. Repeatability of the topographies from phantom measurements: (a) The repeatability of the mean and (b) the repeatability of the aberrations. The repeatability for the moved phantom is given without motion compensation and without outlier removal (blue), without motion correction but with outlier removal (orange) and with motion compensation and outlier removal (green). The repeatability for the static phantom is shown in red. The 95 % confidence interval for the repeatability of the mean is indicated by the gray bars.

Download Full Size | PDF

 

Fig. 10. (a) The compensated motion (red) and the motion applied to the phantom (blue) for a measurement of the moved phantom. The axial motion obtained by prior in vivo measurements (solid blue) was complemented by an unnatural linear part to obtain periodic motion (dotted blue). (b) The difference between the compensated motion and the applied motion.

Download Full Size | PDF

 

Fig. 11. Effect of the motion compensation on the topography for a measurement of the moved phantom. (a) Theoretical axial curvature map. (b) Axial curvature map without motion compensation. (c) Axial curvature map with motion compensation. (d) The mapped differences between the curvature map with motion compensation and the theoretical axial curvature map. The axial curvature maps are color-coded with increments of 0.5 D (Fig. 7(a)). The difference map is color-coded with increments of 0.1 D (e). The corresponding difference measures can be found in Table 1.

Download Full Size | PDF

Table 1. Effect of the motion compensation on the difference to the theoretical topography for the measurement of the moved phantom shown in Fig. 11.

View Table | View all tables in this article

In vivo repeatability Figure 12 shows the in vivo repeatability of our OCT-based topography with and without motion compensation, together with the repeatabilities of the reference topographies. For the mean of the central zone, our method shows inferior repeatability than both reference topographies. For the other zones, the motion compensation enhances the repeatability to a level that is between the repeatabilities of the Atlas and the Pentacam.

 

Fig. 12. Repeatability of the topographies from in vivo measurements: (a) repeatability of the mean and (b) repeatability of the aberrations without motion compensation (blue), with motion compensation (orange), for the Atlas (green) and for the Pentacam (red).

Download Full Size | PDF

In vivo equivalence Figure 13 shows the difference between our OCT-based topography, the Atlas topography and the Pentacam topography for the individual eyes, using the difference mean and standard deviation of the differences described in Section 2.8.3. Figure 13(a) shows the difference mean over the mean of the mean map. The mean map is generated by averaging our topography and the Atlas topography. Figure 13(b) shows the standard deviation of the differences over the standard deviation of the mean map.

 

Fig. 13. Difference mean (a) and standard deviation of the differences (b) between the devices for the different zones. The mean and standard deviation of the difference mean as well as the mean of the standard deviation of the differences are listed in Table 2.

Download Full Size | PDF

Table 2. Mean and the standard deviation (SD) of the zone differences between our method, Atlas and Pentacam. The $p$-value for the mean of the difference mean is calculated using the one sample t-test for the null hypothesis of zero mean difference.

View Table | View all tables in this article

Table 2 summarises the measures of equivalence (see Section 2.8.4) between the topographies. Although the comparison shows some statistically significant differences, we judge the equivalence from a clinical point of view. To judge the clinical equivalence of our OCT-based topography, we use the difference between the Atlas and the Pentacam as primary baseline as they both are clinically approved topographers. As an additional baseline, eyes are commonly classified as ametropic if their absolute refractive error is above 0.5 D [20]. Further, vision corrections are typically available in 0.25 D steps [21]. This indicates that refractive errors below 0.25 D do not have high clinical relevance.

The mean of the difference mean between our method and the reference topographies is consistent with the difference between the reference topographies. The reference devices show a difference in the central zone, which we did not expect on this order of magnitude, but is still below 0.25 D.

The standard deviation of the difference mean is lower between the reference topographies than between our method and both of the reference topographies, especially for the central zone.

The mean of the standard deviation of the differences, the measure for the equivalence of the aberrations, is on a similar level for all zones.

Comparison with McNabb et al. A direct comparison with the results of McNabb et al. [22] is not possible because they did not evaluate the mean curvature and aberrations for the different zones. Instead, they calculated a spherical equivalent power, combining the anterior and posterior corneal curvature. Because we did not evaluate the posterior corneal surface, we are not able to provide a similar measure. However, McNabb et al. reported the repeatability of the mean anterior radius of curvature over a 5-mm-diameter area. We thus additionally determined the mean curvature over a 5-mm-diameter. The resulting SD of the paired differences is 40 µm in radius or 0.23 D in refractive power, assuming a refractive index difference of 0.3375. McNabb et al. reported a repeatability of 26 µm for their DSOCT method, calculated by the population mean of the SD of repeated measurements. Assuming that paired differences are pulled from normal distributions with the same standard deviation of 26 µm, the expected SD of the paired differences is 26 µm $\times \sqrt {2} =$ 37 µm, which is comparable to our repeatability of 40 µm.

4. Discussion

We presented a complete pipeline for OCT-based corneal topography, including novel methods for scanning, motion compensation and reconstruction. We validated the methods by phantom and in vivo measurements. Unlike others [22], we directly compared axial curvature maps instead of derived values such as mean curvatures or axes. The comparison yields separate measures for mean refractive power and aberrations, allowing a more differentiated assessment. The reported results show that our OCT-based topography is comparable with established topographers in terms of repeatability. The equivalence of the aberrations is also at a similar level compared to the reference topographies. For the equivalence of the mean power, we need to differentiate between the mean equivalence and the variations. The mean equivalence is comparable with the reference topographies, whereas the variations in the mean power difference are increased, especially for the central zone.

Nevertheless, we demonstrated that our motion compensation can handle typical eye movements and enhances the repeatability and equivalence of the topography. The results indicate that motion compensation is crucial to achieve a performance that is comparable with established topographers. Other than current methods, we compensate for motion in three dimensions in a single optimisation step. We want to emphasise, that no part of the scan is considered motion-free as done in common methods. Instead, we use a continuous motion model which enables motion compensation with higher temporal resolution and without discretisation artefacts. We therefore presented, to our knowledge, the first continuous motion compensation for OCT measurements of the cornea.

The presented regularisation scheme addresses the stability problems occurring in the curvature space of the Zernike reconstruction. Results show that the regularisation enhances the repeatability of the corneal topography without sacrificing the equivalence with reference topographers.

We want to highlight that the presented motion compensation and reconstruction methods can easily be applied to other scan systems. We further want to emphasise that although we use the motion compensation together with our optimised scan pattern, the method itself is generic and can be applied to other scan patterns with ease. However, there are some limitations that we want to discuss in the next section.

4.1 Limitations and further research

First of all, we optimised and validated the regularisation with measurements of healthy eyes. Thus, the applicability to highly irregular eyes has to be investigated in further studies.

Additionally, the presented scan pattern features a non-uniform scan point distribution, which is suboptimal regarding the stability of Zernike reconstruction [3]. A uniform distribution could be obtained by modulating the speed along the trajectory. However, this would lead to scanner control signals with larger bandwidth, which is why we did not do so in the first place.

Another point to mention is that the optimal chunk size for segmentation and motion compensation, together with the the order of the motion polynomials, are still under investigation. On the one hand, a smaller chunk size for segmentation could reduce the number of removed A-scans during segmentation and outlier detection. Further, the order of the motion polynomials could be reduced, which would potentially increase the stability of the motion compensation. On the other hand, a smaller chunk size potentially reduces the general stability of the Zernike reconstruction used for plausibility check. Increasing the order of the motion polynomials would increase the maximal frequency of the motion compensation but would in return potentially reduce the stability of the motion compensation.

The axial resolution of our system is also low compared to other systems used for corneal topography [6,8], a limitation resulting from the large axial scan range. However, as mentioned above, our methods can be applied to other scan systems with ease. By using a system with a higher axial resolution, the regularisation of the reconstruction could be reduced while maintaining similar repeatability, enhancing the ability to represent irregularities.

Finally, we only compensate for linear motion of the cornea. The compensated motion includes the linear motion that is induced by the head movements as well as the linear corneal motion induced by rotation of the eye. The results show that the compensation of the linear motion is crucial for the resulting topography. However, there is a rotational part of the motion that is not handled by our motion compensation. The investigation of the impact of rotational motion and the inclusion of the rotational motion in our model is the topic of current and future research.