1. Introduction

Precise measurement of the cornea’s refractive power is crucial for the reliable planning of refractive and cataract surgeries, since the cornea accounts for most of the eye’s refractive power. Corneal topography provides two-dimensional maps of the cornea’s curvature, which enables the assessment of its refractive power and aberrations. Therefore, the accuracy of the corneal topography has direct influence on the reliability of diagnoses and surgical outcomes. An accurate determination of the corneal shape allows to detect pathological deformations, such as keratoconus, at an early stage [1]. Corneal topography also enables geometrically correct measurement of sub-surface structures.

Placido disc reflection, scanning slit and Scheimpflug photography are established methods for corneal topography and represent the current clinical standard. Optical Coherence Tomography (OCT) is commonly used for cross-sectional or volumetric imaging, while its use for corneal topography is still limited. While the traditional photography-based principles enable the acquisition of 2D or 3D geometrical information of the cornea at once, OCT relies on the sequential gathering of one-dimensional axial depth profiles (A-scans) to assemble two-dimensional scans (B-scans) or volumetric scans. This increases the measurement duration and hence increases the susceptibility to disturbances like blinking of the eyelid or abrupt eye movements. These disturbances can result in distorted or incomplete scans that require repetition of the measurement. Repetition increases the examination duration and results in stress for the patients. The probability of a distortion-free measurement can be low, especially for patients suffering from tremors.

For corneal topography by OCT, the central area of the cornea is scanned at discrete locations distributed over the area, and at each location the distance to the anterior corneal surface is obtained by image segmentation. Topography maps can then be generated from a parametric representation of the reconstructed surface. Although there is no consensus about the ability of Zernike polynomials to represent all visually significant aberrations [2], they are state of the art for the description of optical aberrations and the reconstruction of optical surfaces. Due to the fact that Zernike polynomials lose their orthogonality in the discrete case [3], the used scan pattern has a direct influence on the problem’s numerical condition. To achieve an accurate Zernike reconstruction, the scan pattern has to account for the radial and tangential oscillations of the Zernike polynomials. With increasing Zernike orders, the radial and tangential frequencies of the Zernike polynomials increase and require finer radial and tangential sampling. Finer sampling improves the numerical condition and potentially enables higher order reconstruction. The practicality and duration of a scan pattern are limited by the A-scan rate and scanner dynamics of the used OCT system.

To date, no scan pattern is reported to feature a balanced scan point distribution, refine over measurement time and enable reliable corneal topography from partial scans. Classic raster scanning is mainly used for volumetric imaging but also for the measurement of the cornea with subsequent Zernike reconstruction [4]. Although complete raster scans usually feature a balanced scan point distribution, missing parts lead to a patchy scan coverage and unreliable reconstruction. Recently, so-called orthogonal scan patterns [5] were used to acquire motion corrected retinal OCT volumes. These scan patterns consist of sequentially acquired raster scans with alternated fast axis. Because the B-scans acquired along the fast axis are almost motion free, the A-scans of the volumes can be registered and realigned retrospectively, resulting in a motion corrected volume. Although some kind of refinement can be achieved by the realignment of the A-scans, this type of scan pattern is not refining per se.

Another commonly used type of scan pattern consists of several meridional B-scans centered around a common point, which is usually aligned to the apex of the cornea. An advantage of these patterns is the simplification of axial motion correction, because all meridional scans share the apex as a common point. A variety of these so called radial scan patterns were used for different geometrical measurements of the cornea [6–8]. McNabb et al. [9] presented a special radial scan pattern for their Distributed Scanning OCT method. The scan pattern consists of 20 radial B-scans with 500 A-scans each. Spatially adjacent A-scans are not gathered sequentially, but distributed over five sub-sampled passages. This method enables finer axial motion correction compared to classic radial patterns. The method was used for the reconstruction of the cornea with 5^th-order Zernike polynomials and has proven its applicability for corneal topography of astigmatism [10]. Although the distributed scan pattern refines over time, this refining is limited to the radial direction. The radial scan patterns, including the distributed scan pattern, feature a concentrated central point density. The resulting unbalanced sampling potentially leads to bad numerical conditioning and reliability of Zernike reconstruction.

In this paper we introduce a new scanning method for OCT based corneal topography. The method features a balanced scan point distribution which refines over measurement time and remains balanced, even when part of the scan is missing. This enables reliable Zernike reconstruction from short or fractional scans and reduces the need for repeating measurements due to distortions i.e. blinking or abrupt movement of the patient.

2. Methods

2.1. Scan pattern

In this paper we introduce a novel scan pattern which is based on Fermat’s spiral and implements two new properties: repetition and golden angle refinement. The use of Fermat’s spiral is advantageous since it combines fast coverage of a circularly shaped area while featuring a balanced scan point distribution. The Fermat spiral was proposed by Navarro et al. [3] for their invertible discrete Zernike transform – providing non-redundant sampling and completeness of the resulting Zernike expansion.

The basic Fermat spiral is defined by the equation r = ±θ^1/2, where r is the radius and θ the angle of the polar coordinates. To implement repetition, the spiral was adapted for a radius that repetitively sweeps linearly from zero to a maximal radius (of typically 3.75 mm) and back, resulting in a spiral that sweeps out and in. One cycle of the resulting pattern consists of two adjacent sweeps: one from the center to the circumference and one back. The first cycle of the scan pattern is shown in Fig. 1(a).

 

Fig. 1 The refining scan pattern with a sampling rate of 2 kHz. (a) Scan points of the first cycle and the initial orientations of the first four cycles. Each successive cycle is rotated with respect to the previous one by an approximation of golden angle. (b) Scan points of the complete pattern consisting of eight cycles.

Download Full Size | PDF

For refinement, a continuous rotation was added to each cycle, such that two subsequent cycles are rotated by an angle θ_r with respect to each other. To obtain a well balanced sampling for any subsection of the scan pattern we used the golden angle θ_g = 0.381966 × 2π as basis for the refinement (see next section). We approximated the golden angle by 0.375 × 2π to obtain a finite scan pattern. As shown in Fig. 1(b), this approximation of the golden angle leads to a repetition after eight cycles, while preserving the advantageous property of balanced refining.

Due to this particular choice of spiral trajectory, the sampling density continuously refines in tangential as well as in radial direction. We defined 256 scan points per cycle by regular sampling in θ which results in 8 × 256 = 2 048 points in total. For a sampling rate of 2 kHz this leads to a duration of 1.024 s.

2.1.1. Golden angle refinement

Our technique uses a continuous rotation of the scan pattern to achieve a refinement of the scanning with time. By using the golden angle as basis, our method uses a property of the golden ratio that is also used by nature. Growing plants position their new leaves with as little overlap with previously grown leaves as possible in order to avoid unwanted shadowing of previous leaves. Some flowers leave an angle between succeeding leaves of around 137.508^° to minimize the overlap [11]. This angle is also known as golden angle and is the smaller angle created by sectioning a circle (360^°) according to the golden ratio. The golden angle in radians is given by

The golden angle is already used for iterative image reconstruction in computed tomography (CT) [12]. To achieve good iterative image reconstruction, each projection has to contain information about the object which is as independent from the information contained by the previously taken projections as possible. Because this optimisation problem is quite similar to the one of the plants, the golden angle was beneficially used as angle shift between subsequent projections.

For our refinement, we used an approximation of the golden angle. Because the golden ratio is an irrational number, subsequent rotation around the golden angle would cause infinite refinement of the scan pattern and would continuously fill the gaps in the angular sampling.

2.2. Evaluation

We evaluate the performance of the method numerically and by measurements with the proposed scanning scheme. To minimise the experimental uncertainties, we rely on glass test surface measurements for the verification of the method. The use of test surfaces is also proposed in technical standards for the assessment of corneal topographers [14]. The practical applicability is shown separately by means of an in vivo use case. The radii of the two used glass surfaces were chosen to match the curvature of two typical corneas with and without astigmatism. A toric surface with radii of 7.585 mm and 7.987 mm was used to represent a cornea with astigmatism. A spherical surface with a radius of 7.805 mm was used to represent an astigmatism-free cornea. The test surfaces were manufactured by Spectros (Ettingen, Switzerland) and kindly provided by Haag-Streit (Köniz, Switzerland). We performed two scans on both of the two glass surfaces. To receive the points for the subsequent reconstruction, the front surface was segmented in the acquired scans.

We employ three measures for the verification of the proposed scanning method. First, the condition number assesses the theoretical bound for the accuracy of the Zernike reconstruction. Second, the standard error of the coefficients assesses the confidence in the reconstructed Zernike polynomials. Furthermore, the root mean square (RMS) of the Zernike coefficients’ difference yields a measure for the difference between two Zernike reconstructions.

For the in vivo use case, we used measurements of a healthy patient acquired over the course of a clinical data acquisition. This data was acquired at the Eye Clinic of the University Hospital Basel with approval of the local ethics committee. The device was aligned to the apex of the eye by using a live view generated from two orthogonal B-scans.

2.2.1. Scan system

We used a custom swept-source OCT system at 1060 nm featuring telecentric scanning and a sweep rate of 30 kHz. A corresponding optical schematic can be found in Fig. 1(b) of the publication of Sarunic et al. [13]. We oversampled the presented scan pattern 15 times to receive a robust segmentation value for the scan pattern points (see next section). This leads to an actual sampling rate of 2 kHz.

The bandwidth-limited full-width at half maximum axial resolution of the system is 40 μm, and the −6 dB distance is larger than 20 mm. Image segmentation is based on data interpolated to 12 μm axial pixel dimension and an A-scan length of 10 mm in air. The power on the sample is 940 μW and the OCT sensitivity is 103 dB.

The scan pattern was calibrated based on scans of a test surface with known physical dimensions. The mean absolute deviation between the calibrated points and the specified points was 37 μm with a maximum of 73 μm. The uncertainty of the calibrated points due to scan-to-scan variations and uncertainties of the calibration method was below 30 μm.

2.2.2. Segmentation

We used a graph-based segmentation method, similar to the one used by LaRocca et al. [15], to derive the points for the following reconstruction. The segmentation was used to identify the position of the front surface in axial direction. To enhance the robustness of the segmentation, the method was adapted to use edges spanning several adjacent A-scans. As a result, one segmentation value was provided for every 15th A-scan. One segmentation value therefore is based on several adjacent A-scans.

2.2.3. Discrete Zernike reconstruction

Given the x- and y-coordinates defined by the scan pattern and the z-coordinates from the segmentation, the cornea front was fitted with Zernike polynomials by linear least squares. The least squares estimation for the highly overdetermined system X β = z was found by solving

2.2.4. Condition number

The condition number is a measure for the sensitivity of the solution to errors in the input. For a linear equation of the form X β = z, the condition number is a property of the problem matrix and is defined by κ(X) = ‖X⁻¹‖ × ‖X‖, or rather by κ(X) = ‖X⁺‖ × ‖X‖ for our overdetermined system, where X⁺ is the Moore-Penrose pseudo inverse of the non-square problem matrix. When using the euclidean norm, the condition number is given by the ratio of minimal and maximal singular values

From this equation one can see that κ(X) ≥ 1. A condition number of one indicates that the error in the input data is not amplified.

2.2.5. Coefficients standard error

The coefficients standard error is strongly related to the confidence bound for the fitted Zernike coefficients [16]. The confidence bound vector is calculated by c = b ± t s where b is the vector with the fitted coefficients and s is the coefficients standard error vector. The factor t can be derived from the t-distribution and the desired confidence level. Because we have a highly overdetermined system, the distribution was assumed to be constant and the coefficients standard error was used as the measure for the confidence. The coefficients standard error vector s is given by the diagonal of the covariance matrix of the coefficients Σ

The estimated covariance matrix of the coefficients is given by

3. Results

3.1. Validation

Figure 2(a) shows the condition number’s dependency on the number of cycles (256 points each) used for Zernike reconstruction. The curves show the course of the condition number for a certain number of points starting at the beginning of the scan. Box plots are produced from the eight sections of the scan with a particular number of adjacent cycles (the last and the first cycles are considered adjacent). As expected, the box plots show very low variance and are therefore only visible as horizontal lines. This means that the numerical condition is only dependant on the number of adjacent cycles and not the specific cycles that are used for reconstruction. The condition of the refining scan pattern steadily decreases and converges to one (blue line) while simple repetition of one cycle does not improve the numerical condition (green line). As a comparison, a radial scan pattern with 20 meridional scans with 500 points each yields a condition number of 2.5. This is worse than the condition of a single cycle of the proposed scan pattern and can not be improved by longer measurement.

 

Fig. 2 (a) Condition number of the Zernike reconstruction depending on the number of adjacent cycles included for the proposed refining pattern (blue) and repeated first cycle (green). (b) Coefficients standard error (blue) and the convergence of the Zernike coefficients (green) for the four test surface reconstructions depending on the number of adjacent cycles included.

Download Full Size | PDF

Analogously, Fig. 2(b) shows the mean coefficients standard error for the four test surface reconstructions and the convergence of the reconstruction result measured by the root mean square of the Zernike coefficients’ difference. The differences are calculated between the Zernike reconstruction, with the points of a given number of adjacent cycles, and the reconstruction of the complete scan. The characteristics are consistent with the course of the numerical condition. The convergence of the reconstruction confirms that cycles can be removed from the measurement with minimal negative impact on the reconstruction reliability. As stated in the introduction, reliable reconstruction is a prerequisite for reliable topography.

3.2. Use case: corneal topography from a partial scan

Figure 3 shows a segmented in vivo scan from a healthy eye with a gross disturbance due to blinking of the patient. The part with the disturbance was manually excluded from the reconstruction. As seen by the small oscillations of the segmented surface, the part after the disturbance is slightly misaligned. To compensate for the misalignment between the parts, we aligned the segmented points of the two parts by means of two separate low order Zernike reconstructions. The shift between the two reconstructions was determined from the Zernike coefficients as described in [17].

 

Fig. 3 The eight individual cycles of a scan with a blinking artifact in the lower row. The segmentation of the anterior cornea is shown in white.

Download Full Size | PDF

Figure 4(a) shows the included scan points and the resulting topography map. After the overall reconstruction with 6^th-order Zernike polynomials, the curvature map was calculated for the central 4 mm disc as described by McNabb et al. [18]. For comparison, Fig. 4(b) shows the curvature map calculated from the reconstruction of a disturbance-free second scan of the same eye. On the used standard scale for corneal topography, both essentially look the same.

 

Fig. 4 Generated corneal topography maps in diopters and the scan points used for Zernike reconstruction. (a) Topography map from a scan with blinking artifact. (b) Topography map from a disturbance free rescan.

Download Full Size | PDF

4. Discussion

The results of the validation show that the scanning method enables well-conditioned reconstruction from fractional scans. As shown, distorted parts of scans can be excluded from reconstruction to enable reliable corneal topography from distorted measurements. With this approach, the need for repeated measurements can be reduced. Because of the repeated coverage of the same area, distortions can easily be detected – either manually or automatically. Although our work is focused on the anterior corneal topography, the advantageous properties also apply to posterior corneal topography or the generation of pachymetry maps.

One could argue, that reconstruction from fractional scans can also be achieved by simple repetition of a basic scan pattern without refinement. However, as we have shown above, repetition does not enhance the numerical condition. Thus without the golden angle refinement, the numerical condition is limited by the basic scan pattern. With the refinement, longer measurement does not only enhance the reliability by increasing the number of points but also by improving the numerical condition. To achieve a similar performance with repetition of a scan pattern, one would have to find a scan pattern that has the same number of points as one cycle of the refining scan pattern while achieving the same numerical condition as several cycles. This is hard to achieve as one cycle of the refining scan pattern is already optimised for good numerical condition itself. Our method therefore eases the trade off between fast coverage and good numerical condition. We want to emphasise that the spiral scanning scheme itself already enables smaller measurement durations than raster scanning and enables higher Zernike orders compared to radial scanning schemes.

The critical reader might pose the question whether subsequent rotation around a fraction of a circle would not boil down to similar capabilities as rotation around the golden angle. This is not the case. Although this approach would also lead to a refinement of the scanning, the conditioning and reliability of partial reconstruction would be limited because of unbalanced sampling. Figure 5(b) illustrates the suboptimal sampling achieved by subsequent rotation around one eighth of a circle scan compared to the golden angle sampling in Fig. 5(a).

 

Fig. 5 The first four cycles from the proposed refining scan pattern (a) with golden angle rotation and an alternative refining scan pattern (b) with an angle shift of $\frac{1}{8}\times 2\pi$ radians.

Download Full Size | PDF

By using a different approximation of the golden angle, the number of cycles, the measurement duration and the overall point density can be adapted. For instance, an increased point density can be required if higher Zernike orders are desired. However, the optimal point density is limited by the expected amplitude of lateral motion.

The main intent of the work was to overcome the corruption of measurements by abrupt disturbances. As we have shown, the presented scanning method yields a solution for this problem. Apart from this, the main characteristics of the scanning method – it enables reconstruction from small fractions of the scan while continuously optimising the overall sampling – also is attractive for other applications. Apart from abrupt disturbances, so called fixational eye motion is another limiting factor. With the presented method, motion could be detected and compensated by cycle-wise or even half-cycle-wise reconstruction and alignment – assuming that motion during one cycle can be neglected.

Disclosures

Joerg Wagner: Haag-Streit Group, Switzerland (F)