Direct modeling of foveal pit morphology from distortion-corrected OCT images

: Inherent distortions aﬀect the spatial geometry of optical coherence tomography (OCT) images and consequently the foveal pit dimensions. Distortion correction provides an accurate anatomical representation of the retinal shape. A novel approach that automatically extracts foveal pit metrics from distortion-corrected OCT images using a sum of Gaussian function is presented. Foveal width, depth and slope were determined in 292 eyes with low ﬁtting errors and high repeatability. Comparisons to undistorted scans revealed signiﬁcant diﬀerences. To conclude, the internal OCT distortions aﬀect the measurements of the foveal pit with their correction providing further insights into the role of foveal morphology in retinal pathologies and refractive development.


Introduction
The development of optical coherence tomography (OCT) has led to cross-sectional in-vivo imaging of the anterior and posterior segment. Given the importance of central vision, characterizing foveal structure in normal and diseased retina is of great interest to research and clinical application. Foveal morphology can be described by the shape of the foveal pit, which is characterized by the histological absence of the inner retinal layers [1]. Altered foveal morphology plays a role in the diagnosis and classification of different pathologies, such as Parkinson disease [2], blue-cone monochromacy [3], albinism [4,5], retinopathy of prematurity [6] and other retinal diseases [7]. Moreover, previously published studies evaluated differences in normal foveal shape based on fixation preferences and cone density [8], refractive error, gender or ethnicity [3,9].
Already in 1980, Williams [10] concluded from psychophysically measured differences in directional sensitivity that the foveal pit follows a simple Gaussian curve. However, modern OCT imaging enabled a more direct and accurate analysis of the foveal shape. The initial measurements of foveal parameters on OCT images were performed manually [11], followed by automated analyses [4,6]. Using a symmetrical Difference of Gaussians function, Dubis et al. [3] firstly described the foveal contour by a mathematical model. Wagner-Schumann et al. [9] utilized the same approach with an additional magnification correction factor for the difference between the nominal axial length given by the OCT device and the measured axial length of the respective eye. Later, other mathematical models were evaluated, which took the asymmetry of the foveal pit -mainly due to the nasal thickening of the nerve fiber layer [12] -into consideration [2,13,14]. Most studies analyzed the macular retinal thickness maps to indirectly describe the outer border of the inner limiting membrane as foveal contour [2,3,7,9,13,14], whilst others extracted manually the parameters like width, depth and slope from cross-sectional B-scans [11]. Nevertheless, only the automated and direct image-based analysis of HD OCT images along with the description of the foveal contour as a mathematical function allows the objective quantification of the fovea.
It must be noted, however, that the resulting B-scan images from the OCT measurements get flattened during the assembly of the individual A-scan lines and the image post-processing, and therefore does not represent the true spatial geometry of the originally curved retina, as already recognized and corrected by previous study groups [15][16][17]. This fact implies that the actual foveal pit morphology should be analyzed from distortion-corrected OCT scans. Therefore, the aim of the current study was to develop a mathematical function that describes the foveal contour of distortion-corrected OCT images and automatically extracts the parameters of interest, like width, depth and slope. The model fit was assessed, validated and compared to a previously described algorithm using the retinal thickness profile from uncorrected OCT scans [3,9].

Study participants
In total, n = 292 eyes from 146 self-reported healthy participants of Caucasian origin with a mean age of 24 ± 3 years were included in the analysis. Participants with ocular pathologies, corneal laser surgery or insufficient OCT signal strength < 6 and thus reduced image quality of the OCT scans were excluded. The mean spherical equivalent refractive errors of the each 146 right and left eyes were −2.19 ± 1.87 D and −2.1 ± 1.78 D. Axial lengths ranged from 20.84 mm to 26.52 mm, with an average of 24.14 ± 1.03 mm and 24.11 ± 1.04 mm. To analyze the repeatability of the modeling of foveal morphology, data of a subset of 10 right eyes was evaluated. The retrospective study adhered to the Declaration of Helsinki and was carried out at the Institute for Ophthalmic Research in Tuebingen.

Instrumentation and measurement procedure
Each participant underwent two different OCT scan patterns with a spectral-domain OCT device (ZEISS CIRRUS HD-OCT 5000, Carl Zeiss Meditec Inc., Dublin, CA, USA), which included a 512x128 volume scan as well as one horizontal HD line scan oriented on the foveal center. Moreover, one measurement of axial length with a swept-source based ocular biometer (ZEISS IOLMaster 700, Carl Zeiss Meditec AG, Jena, Germany) was recorded. Objective refraction was performed using a wavefront aberrometer (ZEISS i.Profiler plus, Carl Zeiss Vision GmbH, Aalen, Germany). The subset of participants for the repeatability analysis underwent the same measurement procedure twice.

Modeling of the foveal morphology
The analysis for the foveal pit morphology followed three main steps: (1) distortion correction of the OCT scan image, (2) edge detection of the inner limiting membrane in the corrected image as the foveal contour and (3) fitting of the mathematical function to the contour with extraction of the foveal pit parameters. Image processing and model fitting was performed with custom-written software in MATLAB (MATLAB 2017b, The MathWorks, Inc., Natick, MA, USA) and OpticStudio (OpticStudio; Zemax, LLC, Kirkland, WA, USA). Distortion correction was achieved using a validated algorithm by Straub & Steidle [17], which estimates the retinal radius of curvature and accordingly corrects the HD line scan. The algorithm performs a three-dimensional ray tracing through the OCT scan optics and ocular surfaces of a customizable optical computer model where individual eye and scan parameters can be modified. In the current analysis, the axial length was adjusted for the individual participants. The retinal shape was estimated by calculating the location of the line of constant group delay across the curved retina. Knowing the numerical displacement of the line, the distortion-corrected OCT image can be reconstructed as shown in Fig. 1.
Afterwards, the corrected images were pre-processed in the following steps in order to detect the inner limiting membrane (ILM) reliably. The scan image was cropped laterally to remove artifacts without image content caused by the distortion correction, which would lead to an erroneous edge detection. Afterwards, a 3x3 median filter was implemented to smooth the image followed by conversion to a binary presentation by applying a grey value threshold level of 50 and a Sobel edge detection algorithm [18]. Using a dilation function with a disk-radius of 2 units, discontinuous edge points were dilated and connected to each other. The number of discontinuous edge points was mainly dependent on the contrast of the B-scan and thus the quality of the detection of the ILM separated to the adjacent retinal layers. An area opening function removed all connected components of the detected edges that had fewer than 1000 pixels. Subsequently, a Sum of Gaussians (SoG) function with three summed terms and each consisting of three fitting parameters a, b, and c was fitted to the central 1000 pixels of the ILM contour (Equation 1).
In order to automatically extract the foveal pit metrics, the squared second derivative (see Fig. 2) from the fitted model was calculated. Local maxima of the resulting function were determined and sorted in a descending manner along with the peak locations. The peak with the highest amplitude represents the foveal center C, while the two successive peaks describe the points A and E on either side of the foveal center with maximum curvature of the ILM, as shown in Fig. 2. The Euclidian distance between A and E was defined as foveal width. The foveal depth was calculated from the perpendicular distance from point C to the connecting line A-E. The maximal foveal slope was extracted via the arctan of the minimum and maximum of the first derivative of Equation 1, representing points B and D on Fig. 2. The temporal and nasal slope were averaged to obtain a single value for the slope. For each of the metric measurements, the isotropic pixels were converted into millimeters by a factor of 2.4 µm/pixel. All scan images were checked for misdetection of the ILM contour and were corrected if needed. Semi-automated correction was performed by either changing the median filter or the threshold value. The subsequent automated curve fitting and extraction of foveal pit metrics were unaffected from image pre-processing changes in any case. Additionally, the results from the developed fitting model using distortion-corrected images were compared to the results from the collected data analyzed by the algorithm from Dubis et al. and Wagner-Schumann et al. [3,9]. Their published algorithm, using a Difference of Gaussian (DoG) fit, was applied to the Macula Cube scan from the current dataset. Therefore, the central horizontal profile was extracted from the retinal thickness map, which uses the relative distance between RPE and ILM, and the sampling rate was adjusted to match the line scan underlying the SoG fitting.

Data analysis
For the assessment of the intrasession repeatability for the foveal width, depth and slope, the Intraclass Correlation Coefficients ICC(2, k) and coefficient of repeatability were calculated out of two consecutive measurements [19]. The coefficient of repeatability is described as the standard deviation of the differences between the test-retest measurements, multiplied by factor 1.96 [20]. Bland-Altman analysis and the Wilcoxon signed rank test were additionally performed to compare and evaluate the agreement of the foveal pit metrics extracted from distorted versus undistorted images from the current and previous published algorithm [3,9].

Assessing model fit accuracy
The average root mean square error (RMSE) between the fitted SoG function and the detected ILM is as low as 2.16 ± 1.94 µm and 1.89 ± 1.35 µm for the right and left eye. The repeatability analysis showed an excellent correlation with ICCs > 0.99 for the morphological parameters, as shown in Table 1. The coefficients of repeatability were 35.47 µm, 3.75 µm and 0.20 • for the width, depth and slope analysis, respectively. Corrections of the ILM contour detection were necessary in 16.94 % of the images.

Foveal pit parameters
The extracted foveal width, depth and slope with associated repeatability metrics are shown in Table 1. All foveal pit parameters are characterized by high intersubject variability as indicated by the standard deviation and total range. Pearson correlation showed a strong correlation between both eyes for all analyzed parameters.

Comparison to foveal pit morphology from distorted images
A total of 113 right eyes could be evaluated with both algorithms. The maximum foveal slope is the only common and thus directly comparable parameter between the newly developed algorithm with distortion-corrected scan images and the previous algorithm using distorted scans [3,9]. The Wilcoxon signed rank test revealed a significant difference (p < 0.001) between the results of the fitted two models -DoG with uncorrected retinal thickness profiles vs. SoG with corrected images. Further Bland-Altman analysis showed a mean difference of 4.17 • between both algorithms and 95 % limits of agreement of 1.74 • and 6.61 • . As shown in Fig. 3, the bias between distorted and undistorted results increased with steeper mean slopes. A comparison between foveal width and depth between the distortion-corrected SoG and the uncorrected DoG fit is possible, however, has to be taken with caution due to the different mathematical definitions of the function landmarks from which the foveal pit metrics were extracted. Foveal pit widths were significantly different from each other with 2043.83±355.97 µm and 1195.13 ± 223.90 µm (p < 0.001) for the DoG versus SoG fit. In contrast, the calculation of foveal depth revealed no difference (p = 0.37) between both methods with values of 107.51 ± 23.48 µm (uncorrected) and 112.50 ± 24.540 µm (corrected). Moreover, the presented SoG fit to the distortion-corrected OCT scans showed a significantly lower average RMSE (2.18 ± 2.00 µm) than the DoG fit to the retinal thickness profile (8.79 ± 2.61µm) on a level of p < 0.001.

Discussion
The developed model for the description of foveal pit morphology consisted of two automated steps, which are firstly the distortion-correction of B-scans and secondly the fitting of a SoG function to the distortion-corrected foveal contour with subsequent extraction of the foveal pit metrics from mathematically defined locations. The utilized SoG function resulted in a low RMSE and thus is able to describe the foveal morphology closer to anatomical proportions as represented by the distortion-corrected images than a DoG function [3,9]. Moreover, the current study found a high intersubject variability in all parameters. To the best knowledge, this is the first automated and image-based evaluation algorithm for foveal morphology from distortion-corrected OCT images.

Advantages of foveal pit morphology modeling from distortion-corrected OCT images
This approach is advantageous in multiple ways.
(1) The corrected B-scans reflect the retinal size proportions closer to the actual anatomical dimensions. Moreover, the algorithm by Steidle & Straub was previously validated using the same OCT on an artifical eye whose parameters were exactly known. They found no deviation between the estimated retinal shape by their algorithm and the true retinal shape of the test eye [17]. In contrast, regular post-processing from the OCT system results in a flattened, thus distorted, representation of the scan area, which subsequently affects the absolute measures of foveal pit metrics [15][16][17].
(2) The correction of the inherent distortion considers the axial length of the respective eye.
The influence of axial length on absolute retinal scan dimensions had not always been considered in previously conducted studies. Either no correction factor was applied [3,6,7,11,13,14] or participants with high hyperopia or myopia and thus short or long axial lengths were excluded [2]. However, as already shown the axial length significantly influences foveal width and slope measurements, but not necessarily the foveal depth [9]. Therefore, individual axial length measurements were applied already during the distortion correction in the current study. In order to increase the precision beyond the measurements of axial length, the distortion correction model considers all optical elements from the OCT system, the individual's measured axial length and assumed values for remaining eye variables using the Arizona Eye Model [17]. This approach shows greater benefit in subjects with extreme axial lengths compared to a more simplified correction factor for the magnification of the scan.
(3) The foveal pit modeling can be performed directly on the exported B-scan images. This approach does not depend on the retinal thickness profiles, as it was the case in most previous studies [2,3,6,9,13,14]. Subsequently, the direct modeling does not only eliminate one additional step but also an additional source of errors, in case of erroneous calculation of the retinal thickness. This could be the case in eyes with retinal diseases with misdetection of the ILM and retinal pigment epithelium for the generation of the retinal thickness maps [21]. Moreover, the purely image-based analysis automatically extracts automatically the points A to E were defined via landmarks in the SoG model and its derivatives. Therefore, inter-examiner influences and manual measurement errors do not need to be considered for this study.
(4) The broad variety in foveal morphology, as well as the e.g. naso-temporal pit asymmetry [12] cannot be described by symmetrical functions alone. The SoG model with three Gaussian terms is mathematically able to describe asymmetric pit shapes, as well as foveae with flat or steep bottoms. In contrast, the DoG function does not reliably fit the undistorted foveal shape, as its derivative does not completely comprise three roots for the points A, C and E. Moreover, the current study reported a RMSE, the difference between the fit result and the ILM contour, of 2 µm, thus, the curve fitting can be neglected as potential source of error. Other studies using asymmetrical algorithms also reported small RMSEs between 1.40 µm [13] and 4.25 µm [14], whereas symmetrical models showed higher deviations up to on average 14 µm [3,9].

Comparison with previous literature and undistorted OCT scans
However, it is difficult to compare the current morphology findings regarding to past studies, due to the great variability of applied analysis methods, mathematical models and definitions of foveal landmarks. Foveal width metrics appear to be mostly affected by the methodological inconsistency across the literature with average values ranging approximately from 1200 µm up to 2500 µm. For example, Nolan et al. [11] defined foveal width in two different ways for their manual measurements, from the point where the nerve fiber layer begins and additionally from the highest crest points on both sides of the foveal pit. These findings are most identical to the current results with 1200 µm and 1400 µm for both definitions, respectively. Automated analysis generally led to a broader foveal width measurements from 1800 µm to 2000 µm with a DoG function [3,9], 2100 µm with a sloped piecemeal Gaussian model [14] and 2400 µm from the points of maximum macular thickness [7]. However, even when the measurement points are similarly defined, such as Nolan et al. [11] and Tick et al. [7], who both measured from the foveal peaks on each side, still found significantly different results for the foveal width. This might be attributed to the use of different OCT types (time-domain vs. spectral-domain) with consequently different scan quality, optically induced distortions, different measurement methods (manually vs. automated) or the image analysis modality (exported vs. within the device software). In contrast, foveal depth and slope results are more consistent despite different methodology, ranging from 110 − 130 µmm and 10 − 12 • . The newly developed SoG model revealed similar results for the foveal depth, whereas the slope values were significantly steeper with an average of 14 ± 2 • . The results for the foveal slopes from the DoG and SoG model were compared to further investigate the influence of the analysis of distorted vs. undistorted images. Using distortioncorrected scans resulted in significantly steeper slopes. Moreover, the bias became more prominent with increasing foveal sloping. This difference can be explained with the general curvature of the distortion-corrected OCT, which affects the slope proportions. Width and depth were not directly comparable between the two algorithms. Nevertheless, it can be expected that the generally small foveal width in the current study also partly results from the "curving/bending" of the undistorted scan image. The use of two different mathematical models (SoG vs. DoG), and the additional consideration of the OCT and eye optics, however, might have been attributed to these differences as well.

Limitations
Ground truth foveal pit metrics in the form of measured parameters derived from histological measurements on extracted human eyes do not exist, mainly due to the difficulties and distortions during tissue extraction and preparation. Using a semi-automated image-based software for histological quantification of foveal pit metrics is only possible for histological preparations from animals so far [22]. However, the aforementioned validation study of the used algorithm [17] on an artificial test eye with known ocular parameters holds strong proof for the accuracy of the used distortion correction. Moreover, the conducted study is bound to the OCT device mentioned in the methods section (ZEISS CIRRUS HD-OCT 5000, Carl Zeiss Meditec Inc., Dublin, CA, USA) as the algorithm only considers the optical components from this device. However, it was extended more currently to a swept-source OCT (PLEXElite 9000, Carl Zeiss Meditec Inc., Dublin, CA, USA). The results from the current study show that the usage of ray-tracing based optical distortion correction for OCT B-scans have clinical implication and thus could be implemented into other types of OCTs for the inclusion into daily ophthalmological and optometric routine.
Another possible limitation of the conducted study is the mainly myopic set of participants. However, a previous study did not report any significant differences in foveal pit morphology between myopes and emmetropes [3]. Moreover, semi-automated corrections might be necessary in case of an erroneous ILM detection, which can occur for example in the presence of prominent vitreal floaters and ILM pathologies, such as epiretinal membranes [23]. The current study performed and evaluated only horizontal scan images, however, there is evidence that the fovea has different horizontal and vertical dimensions [14], which have not been considered yet.
Another drawback lies in the control of the various variables of the optical simulation model, such as the working distance between the eye and the OCT. The working distance influences the optical eye model with ray-tracing and subsequently the estimation of the retinal shape underlying the final correction of the image. The same pitfall accounts for individual biometrical variables, such as anterior chamber depth, gradient refractive index of the lens or corneal curvature. However, the present study did not evaluate the extent between the error in the retinal radius calculation and the corrected image proportions, which were eventually used for the foveal pit metrics. It is assumed that they have a greater impact on the peripheral sections of the OCT image rather than the central part, which was used for the SoG fitting. Moreover, the calculation of the slope inherits an additional possible source of error, as it is obtained via the arctan of the mathematical slope. This requires a reference coordinate system, which was the horizontal x-axis of the image in the current study. However, a tilt of the scan image then leads to an overor underestimation of the slope.

Conclusion
The study provided evidence that a SoG function is able to precisely describe foveal pit morphology in distortion-corrected OCT scans and provided measures of foveal parameters closer to retinal dimensions. It further showed that the correction of distortions has an effect on foveal width and maximal slope. Foveal morphology is of great interest in the screening, diagnosis or function testing of healthy and pathologic eyes. A physiological relationship between foveal morphology and the Stiles-Crawford effect [24] has already been suggested previously [10] but has been re-adressed recently again in the field of myopia research [25]. In this context, future studies could further investigate the foveal pit metrics in different scan directions, such as vertical and oblique angles to analyze ocular growth models in children [26]. Moreover, correlations of these parameters with biometrical and refractive error data could give important insight into foveal morphology associated with axial elongation, together with subsequent transfer of the gained knowledge into clinical practice, such as an additional factor for the estimation of myopia progression and for the efficacy control of the myopia treatment strategies.