Ultrawide field, distortion-corrected ocular shape estimation with MHz optical coherence tomography (OCT)

: Ocular deformation may be associated with biomechanical alterations in the structures of the eye, especially the cornea and sclera in conditions such as keratoconus, congenital glaucoma, and pathological myopia. Here, we propose a method to estimate ocular shape using an ultrawide field MHz swept-source optical coherence tomography (SS-OCT) with a Fourier Domain Mode-Locked (FDML) laser and distortion correction of the images. The ocular biometrics for distortion correction was collected by an IOLMaster 700, and localized Gaussian curvature was proposed to quantify the ocular curvature covering a field-of-view up to 65° × 62°. We achieved repeatable curvature shape measurements (intraclass coefficient = 0.88 ± 0.06) and demonstrated its applicability in a pilot study with individuals (N = 11) with various degrees of myopia.


Introduction
Ocular deformation may be associated with biomechanical alterations in the structures of the eye, especially the cornea and sclera in conditions such as keratoconus [1,2] and pathologic myopia [3,4]. Severe alterations in ocular shape can lead to pathological changes in ocular structures, for example, the formation of posterior staphyloma in pathological myopia, where the outpouching of the sclera leads to abnormalities of the retinal pigment epithelium (RPE), retina, and choroid. Thus, monitoring the changes in the ocular shape of myopic eyes may be useful in detecting early structural changes, before damage and visual impairment occur with the development of actual staphyloma.
However, imaging the entire eye to assess its ocular shape is currently problematic [4]. Conventionally, magnetic resonance imaging (MRI) may be used [5][6][7][8][9], but this approach has many disadvantages: it is bulky, costly, time-consuming, yields low, anisotropic resolution, and uncomfortable and/or contraindicated for some patients. As MRI machines are not readily accessible in many eye clinics, this modality is not currently used as part of routine ocular imaging and assessment. Other approaches, such as ultrasonography [10,11], partial coherence interferometer [12,13], and peripheral refraction [14], have limited sampling rates or resolutions. Optical coherence tomography (OCT) can provide a volumetric representation of the posterior pole with optical resolution and an acquisition time within seconds [15]. In particular, swept-source OCT (SS-OCT) provides a better signal roll-off and larger imaging depth than spectral-domain OCT. With the recent development of megahertz (MHz) swept-source engines, including Fourier domain mode-locked (FMDL) lasers and vertical-cavity surface-emitting lasers (VCSEL), a dense-sampling volumetric scan can be acquired within a cardiac cycle.
A commercial ultrawide field SS-OCT system (Canon, Japan) can scan a field of view up to 23 × 20 mm 2 on the retina, thus it enables detection of wide staphyloma [16][17][18]. However, it remains challenging to measure the actual changes of the ocular shape because of the optical distortions that are inherent to OCT imaging: several measures, such as the distance between machine and eye for a given exam can significantly alter the ocular shape displayed in OCT images. To solve this problem, McNabb et al. developed a whole eye system to simultaneously image the anterior and posterior eye [19]. A ray-tracing model was then used to reconstruct the true retinal shape. Later, the same group proposed a sparse radial scan pattern to cover a disk area 12 mm in diameter at the posterior pole. However, limited by the system speed (100/200kHz), the curvature can be quantified only along each radial B-scan, and a tomogram of the retinal curvature was generated via interpolating the radial scans [20]. Additionally, a similar approach was also applied to commercial devices and yielded comparable results to MRI scans [21], but the knowledge of the optical parameters was not fully known, making inter-scan or inter-patient comparisons difficult.
Based on these limitations, we propose a method to estimate the ocular curvature using a MHz SS-OCT with an FDML laser and an ultrawide field optical design. The ocular biometrics for distortion correction was collected by IOL Master 700, and a localized Gaussian curvature was proposed to quantify the ocular curvature up to 65°×62°. We also tested its repeatability and demonstrated its applicability in a pilot study with individuals (N = 11) with various degrees of myopia. The ability to quantify ocular shape can allow for a more objective staging of posterior staphyloma in pathological myopia and its progression.

OCT system and imaging protocols
An optical schematic of the MHz FDML SS-OCT system (Optores, Munich, Germany) is shown in Fig. 1, and a photograph of the patient interface module is shown in the inset. Briefly, the light source is an FDML laser (central wavelength: 1061 nm, bandwidth: 75 nm; Optores, Munich, Germany) operated at a 1.66 MHz A-line rate. A calibration arm was integrated to generate a monotone with the reference arm for the linear-k resampling. A Hilbert transform was applied to the interference fringes, where the nonlinearity of the phase was resampled by linear interpolation. The switch between calibration and eye imaging was controlled by an electrical shutter in the calibration arm. The coefficients of the k-space resampling were applied to the retinal images, and the coefficients remained relatively stable throughout the imaging session. The sample arm consisted of a collimator, a pair of galvo scanners (6215H, Cambridge Technology, MA, USA), and customized telescope lenses (f 1 = 91.0 mm, f 2 = 42.0 mm) intended to minimize the eye aberration and light transmission at a large oblique incident angle [21]. The ocular lens design catered for ultrawide imaging up to 100°in emmetropic eyes, but a smaller FOV is more suitable for myopia. An iris camera and a fixation target (λ = 635 nm, P = 7µW) were integrated for eye alignment and fixation, respectively. A dual-channel balanced photodetector (WL-BPD1GA, Wieserlabs, Munich, Germany) and a PCIe based digitizer (ATS9373, Alazar Technologies Inc., Montreal, Canada) were employed for fringe detection and digitization. The digitizer operated at 3.68 GS/s clocked internally, sampling 2304 points for each sweep cycle, providing an imaging range of 4.2 mm and >6 mm 6-dB signal roll-off. The sensitivity roll-off was dominated by the limited bandwidth of the photodector [22]. The axial resolution was 4.7 µm in tissue, and the lateral resolution was approximately 45 µm measured on a USAF resolution target (Thorlabs, Newton, NJ, USA) in a model eye (E-106-1, Epipole, UK). The beam size at the pupil was 1.3 mm. The signal-to-noise ratio was 95 dB with an incident power of 3.0 mW onto the cornea, which together with the power from the fixation target was below the maximum permissible exposure (ANSI Z80. . A data acquisition software (Optores, Germany) with real-time B-scan display was used for in-vivo image alignment. For the scan protocol, the field of view (FOV) was 65°in the fast scan direction and 62°in the slow scan direction. The calibration of the FOV was conducted based on a customized phantom eye with a 3D printed bowl-shaped plastic target (next section). Each volume consisted of 1536 A-scans and 993 B-scans. The acquisition time for each volume was approximately 1.6s, with a duty cycle of 60% on the fast scanner. The image was oriented such that deeper ocular structures were closer to DC to improve sensitivity in the macular region.
A total of 11 eyes from 11 normal subjects (34.8 ± 8.7 years old, range: 25-56 years old) with different degrees of myopia were recruited from Singapore National Eye Center and imaged in triplicate without mydriasis and with a 1-minute interval between images. The room light was turned on to constrict the pupil size for better alignment of the pivotal point with respect to the pupil. Ocular biometrics, including axial eye length, central corneal thickness, and corneal curvature were later acquired by IOLMaster 700 (Zeiss Meditec, CA, USA).

System calibration
A customized eye phantom ( Fig. 2) was utilized to calibrate the scan FOV. It consisted of a 3D printed iris (3 mm diameter), an achromatic doublet (AC127-025-C, Thorlabs, USA), and a plastic target with a bow-shape surface, 3D printed from a stereolithography printer (Form 3, FormLab, USA). The curvature of the surface was designed to be 13.00 mm, with four rows of markers extruding 100 µm vertically from the surface. The interval between the markers was 1 mm in the transverse direction. A spacer was placed in between the doublet and the plastic target to place the target at the focal plane of the lens. The plastic surface was imaged with the beam pivot carefully placed at the iris center. The raster scan directions were aligned to two orthogonal rows of markers. The image volume was then compressed onto an elevation map by identifying the location of the maximal intensity along each A-scan. Loci of the markers were pinned manually from the map, and their corresponding scanning angles θ n were calculated as tan −1 (X n /Y n ), where X n and Y n were the horizontal and vertical distances from the nth marker to the pivot, respectively. Subsequently, the relationship between the scanning angle and marker loci showed that the scan location was linearly correlated with the scan angle and minimally affected by the galvo scanner jitter. By extrapolating the linear relation to the image boundaries, the FOV of the scanning protocol was determined.

Distortion correction via ray tracing
The distortion was corrected by modeling the beam propagation in the sample arm (OpticStudio, ZEMAX). The detailed method has been described elsewhere [19,23]. Briefly, optical parameters, including the surface, material, and distance information from galvo scanners, telescope lenses, and human eyes were incorporated into the model, and a Polans eye model [24] was utilized to represent the human eye, where corneal thickness, anterior chamber depth and axial eye length were measured from the IOLMaster 700. The ray-tracing model allowed for simulating the beam scanning by rotating the x-and y-galvo scanners. The optical path length (OPL) between the x-galvo scanner and the retinal surface was then calculated for each incident angle. A total of 300 × 200 incident angles along with the fast and slow directions respectively were stimulated for each eye.
The models of one short eyeball and one long eyeball are shown in Fig. 3. Blue lines with arrows indicate the chief ray propagation to the retina, with one normal to the eye and two oblique ones representing the boundaries of the scanning FOV. It is apparent that identical scanning angles yielded a larger scanning area on a longer eyeball.

Image post-processing
The RPE layer over the entire volume was segmented automatically based on graph theory and dynamic programming [25,26]. The RPE layer on a central scan was initially segmented, and a band 20 pixels above and underneath the RPE layer was utilized as a region of interest for the neighboring scans. This procedure was iterated and propagated to the entire volume. Although the MHz rate of the SS-OCT engine effectively reduces the volumetric acquisition time, motion artifacts were still present that can affect high-order ocular shape descriptors. These motion artifacts, originated from a variety of linear and non-linear motions, including translation, rotation, and pulsatile expansion/contraction were difficult to isolate. To address these artifacts, we propose a simplified motion estimation and compensation method. First, before distortion correction, profiles of the RPE elevation along the slow scan direction were extracted, at the locations 3 mm nasally and temporally eccentric to the fovea. Second, a 7 th order polynomial function fits each profile with a bisquare optimization, where the residuals were calculated (L 1 and L 2 in Fig. 4(A)). Last, a residual map was generated by linearly fitting L 1 &L 2 on each slow-scan position (Fig. 4(A)). The black arrow and red arrow labeled the artifacts that most likely represent a translation and a tilt of the eyeball, respectively. The motion artifacts (residues) were removed from the surface, the result of which was then to correct the distortion.
To compensate for any remaining segmentation inaccuracy, we further low pass filtered the retinal surface with 2D polynomial fitting to the 9 th order, with bisquare optimization. As the ocular shape change from pathology is mostly on a large scale, the high order polynomial fitting only removed the small-scale surface variation while preserving the localized curvature change. The fitting residual was mostly within 20µm and partially came from the inaccurate segmentation underneath the large blood vessels in the superficial plexus ( Fig. 4(B)).

Curvature quantification
To quantify local retinal shape, we adopted the Gaussian curvature (K), which is an invariant topological parameter previously used to describe the stretch or compression of biological tissue or cells [27]. The Gaussian curvature at point p is calculated as the production of the principal curvatures k 1 and k 2 from two normal sections. An illustration of Gaussian curvature calculation on a saddle function (Z = X 2 − Y 2 ) is shown in Fig. 4(C). The Gaussian function K at the saddle point is less than 0 with k 1 >0 and k 2 <0.
The effects of motion compensation and 2D polynomial fitting on the Gaussian curvature are shown in Fig. 4(D)-(F). Without any compensation, the Gaussian curvature map was deteriorated by motion artifacts, appearing as horizontal stripes (Black arrow, Fig. 4(D)) across the entire FOV. After motion compensation, the horizontal lines were largely removed, but curvature variation from the residual motion artifact or segmentation inaccuracies still caused erroneous curvature detection. Notably, some of the errors followed the pattern of the large vessels from the superficial plexus(white arrow, Fig. 4(E)), and the high variation in the fovea was because of a segmentation error (Fig. 4(E)). Further fitting the surface removed local variation while preserving the high-order ocular shape (Fig. 4(F)). A customized colormap [20] was used to enhance the contrast: normal curvature ranges were shown in green, flatter retina was shown in purple, and more concave retina was shown in yellow to orange. A distinct color gradient from purple to pink indicated the transition from concavity to convexity. The entire FOV was separated into five sectors: central 30°sector and four peripheral quadrants (Nasal, Temporal, Superior, Inferior) (Inset in Fig. 4(G)), and the Gaussian curvature was calculated over the entire FOV, as well as in each sector.

Result
Figure 5(A) shows an OCT cross-section of the plastic target. Twenty-one markers from one row could be visualized, and printing errors were present on the right side of the image (white arrow). A volumetric OCT view could visualize all the markers (Fig. 5(B)), which were used for FOV calibration. The distortion-corrected target surface is shown in Fig. 5(C). The radius of the curvature of the best-fitted sphere was 12.87 mm, and the Gaussian curvature was 0.0057 ± 0.00057 mm −2 from all sectors. The designed spherical radius of curvature was 13 mm in SolidWorks, equivalent to Gaussian curvature 0.0059 mm −2 . The detected curvature was well-aligned with the design, and the 2% difference could originate from the limited printing resolution and resin shrinkage during washing and curing.  Figure 6(A) shows a raw 3D rendered posterior pole centered at fovea from a subject (age: 32 years old, Male; Visualization 1). The low-intensity regions in the volume edge (white arrow) were due to the light obstructed by the eyelashes. Figure 6(B)-(C) are cross-sectional scans across the optic nerve head and 15°superior to the fovea, respectively. The images were acquired with the choroid placed close to the zero-delay and flipped vertically for visualization purposes. Red lines indicated the results of RPE segmentation, and an animation flying through the entire volume with the RPE segmentation is shown in Visualization 2.
To test the repeatability of the method, a total of 11 eyes from 11 subjects were imaged three times. Figure 7 shows the representative results from one eye (age = 34 years old, female). Visually, the distortion corrected RPE elevation maps and the Gaussian curvature maps showed clear and repeatable patterns: the ocular surface was flatter (smaller K value) superior to the optic nerve head, and more concave (larger K value) along an arc in the inferior and temporal peripheral regions. Among all eyes, the radius of ocular curvature with the best fitted sphere had perfect repeatability (intraclass coefficient (ICC) = 0.96) and Gaussian curvature in different sectors had good to perfect repeatability (ICC: 0.88 ± 0.06), where the nasal region had the lowest ICC value (0.78) and the central 30°the highest (0.94).
Representative Gaussian curvature maps from four eyes are shown in Fig. 8. No ocular disease was diagnosed among these eyes, and they were arranged with increasing axial length and decreasing spherical equivalent. In subjects 1&2, the Gaussian curvature over the entire FOV was around 0.01 mm −2 , with slightly more concave regions at some of the peripheral quadrants. In subject 4, the retina was flatter in the macular area, shown in purple color with lower Gaussian curvature. The Gaussian curvature from 11 healthy eyes with different spherical equivalent (-1.00D to -8.00D) and axial eye length (23.18 mm to 26.88 mm) were summarized in Table 1. The radius of ocular curvature was 11.26 ± 0.31 mm with the best sphere fitting, equivalent to a Gaussian curvature of 0.0079 mm −2 , and the measured Gaussian curvature was 0.0078 ± 0.0018 mm −2 , equivalent to 11.30 mm radius of ocular curvature. It indicated that Gaussian curvature contains a wealth of information regarding not only the radius of the ocular curvature but also, more importantly, higher-order localized curvature.

Discussion
We demonstrated distortion-corrected ultrawide ocular curvature estimation using a MHz SS-OCT system. A MHz OCT engine could substantially reduce the volumetric acquisition time and subsequently the motion artifacts. A ray-tracing algorithm reconstructed the true ocular shape from OCT scans, and a Gaussian curvature map was proposed here to quantify the localized ocular curvature, which could help us detect macular and peripheral ocular shape abnormalities in various ocular diseases, such as staphyloma in high myopia. Distortion-corrected ocular curvature detection can identify staphyloma more readily. The commercially available OCTs provide a limited FOV that could erroneously call a steep retina contour (within a 6 or 12 mm scan) a staphyloma [28,29]. The machine-to-eye distance could substantially affect the curvature displayed in OCT images: placing the eye away from the pivot could make the retina more curved. Although full knowledge of the optical parameters in the sample arm is required to reconstruct the true ocular curvature, this technique could also detect relative and longitudinal changes of the retinal curvature on standard OCTs if machine-to-eye distance and ocular biometrics are not available.
Full eye imaging modalities, including B scan ultrasonography [10,11] and MRI [5][6][7]9,30] are not commonly used in clinics because of low resolution and high demand for operational expertise. Moreover, MRI is costly, time-consuming, and uncomfortable and/or contraindicated for some patients. In more selective cases where there is already the presence of myopic traction maculopathy, the ultrawide field enabled by the MHz OCT (up to 65°×62°) is also beneficial as it has been postulated that peripheral pathologic traction may sometimes be the driving force that results in pathology that manifests in the posterior pole [31][32][33][34].
The existing hypothesis is that curvature changes in highly myopic eyes may precede the development of outright staphyloma outpouching, and that asymmetry of curvature may be an early indicator of localized weakness that will eventually result in a staphyloma [3]. The Gaussian curvature change is a resultant of the localized biomechanical properties alteration that the value deviates from the normal distribution is a sign of tissue stretching and compressing, instead of bending (Gauss's remarkable theorem) [27]. Some biomarkers could be extracted to discern which highly myopic eyes will remain stable versus which highly myopic eyes will progress to pathology with staphyloma formation, myopic traction maculopathy, or myopic macular degeneration. Previous studies have reported a potential role of choroidal thickness and scleral shape in the development of pathologic myopia and myopia-associated optic neuropathy. Integrating the curvature information with these other metrics, such as choroid thickness, sclera thickness, and their biomechanical properties could help us estimate tissue stress to strain relationship and its role in pathological myopia development in the future [35][36][37].
One limitation of this study is the lack of a precise machine-to-eye distance measure, which as aforementioned, a key parameter in the ray tracing algorithm. Although ultrawide scanning mode without mydriasis allows us to align pivot more accurately (because clipping would occur if the pivot were outside the constricted pupil), the assumption that the pupil center was placed at the pivot could be compromised. The uncertainty in pivotal point alignment could cause 8% and 2% variations in Rc and Gaussian curvature estimations. Potential solutions are incorporating a dual channel for the anterior segment, using an iris camera with a tunable lens to calibrate the machine-to-eye distance, or calibrating the reference arm position during each acquisition. Moreover, while we did not use relay lenses between two scanners, our ray tracing algorithm simulated this non-perfect pupil plane conjugation and corrected for the resulting distortion. Lastly, the vignetting-caused low intensity at the peripheral regions increased the uncertainty of RPE segmentation, where the segmentation errors were manually corrected or excluded for the following quantification.
In conclusion, we propose a method to estimate the ocular curvature using a MHz SS-OCT with a FDML laser and an ultrawide field optical design. High repeatability of the radius of the curvature and Gaussian curvature was achieved. A pilot study tested its applicability on emmetropic and myopic eyes. Future studies will explore its utility in larger myopia cohorts and investigate the ocular biomechanical properties in high myopia.