2.3. Image Processing Pipeline

Using the built-in software of the scanner (Heidelberg Eye Explorer 1.9.10.0, HRA Spectralis Viewing Module 6.16.0) both the inner limiting membrane (ILM) and the Bruch’s membrane (BM) were segmented from each B-scan (Figure 2a). All OCT images as well as their automatic segmentation were visually reviewed and segmentation errors inside a 3 mm radius region were manually corrected upon consensus between two OCT experts (A.M.-G. and I.G.). Lateral scaling, influenced by axial length differences, was automatically adjusted for each subject by the built-in software. As described in [39], the estimation uses the Gullstrand schematic eye model [40] as reference and adjusts the lateral scale based on each subject’s refractive error (considered when the eye is focused during acquisition) and keratometry values. Images and layer segmentation data were exported to vol format for posterior analysis.

Main processing steps of the foveal pit morphology analysis pipeline. (a) Segmentation of the inner limiting membrane (ILM) and the Bruch membrane (BM), (b) total retinal thickness (TRT) calculation, (c) location of the foveal center, (d) fitting mathematical models, (e) computation of geometrical parameters.

All subsequent data processing was carried out using custom software developed in MATLAB 2020b (MathWorks, Inc., Natick, MA, USA). Several external helping functions were used throughout the analysis [41,42,43]. In the first place, the vol files were opened with the OCT Layer Segmentation package of AUtomated Retinal Analysis (AURA) tools [44], an open source MATLAB library for retinal image processing developed in [45]. Coordinates of each A-Scan were retrieved from the vol data and transformed so that the x and y axes represent the temporal to nasal, and inferior to superior directions, respectively. Left eyes were flipped to match right eyes. From the retinal layer segmentation, which included the point-to-point distances from the bottom of each B-scan image to the boundary of ILM and BM, TRT was calculated as:

This step, which is equivalent to performing a flattening of the image where the BM is set as a reference, is helpful to disregard the effect of the eye curvature and to set a common flat reference to compute morphological parameters (Figure 2b). Combining the TRT values obtained for each point (A-scan) of each slice (B-scan) a 2D TRT raw map of the entire surface was obtained for each eye. These TRT maps were used to automatically determine the foveal center and align the scans (Figure 2c). To this aim, four different strategies were implemented and compared (see Section 2.3.1). The location of the foveal center was used to center the TRT maps using a 2D translation.

Finally, centered TRT maps were resampled to two different patterns using triangulation-based 2D cubic interpolation as implemented by griddata function in MATLAB 2020b:

Regular grid of 3 × 3 mm² and a spacing of 0.02 mm. This was used for foveal center location method comparison (see Section 2.4.1).

Radial pattern with 2 mm radius, 24 angular directions and a spacing of 0.02 mm. This was used for morphology analysis and mathematical model comparison (see Section 2.4.2). This was calculated after using only the smooth + min method to locate the foveal center, as it was the method that provided the best alignment.

Using radial data, the pit morphology models, and smoothing methods described in Section 2.3.2 were used to adjust the TRT curves (Figure 2d). Models were adjusted based on the non-linear least squares method as implemented in MATLAB with a maximum number of 1000 iterations and a tolerance of 10⁻⁶ for both the residuals and the model coefficients. The initial values of the coefficients were manually fine-tuned and the option achieving the best results in terms of fitting error was finally used.

Using TRT raw curves as well as the TRT curves obtained after applying the aforementioned approaches, the geometrical parameters described in Figure 2e were computed as follows:

Central foveal thickness (CFT): the TRT value at the foveal center.

Rim height: the point of maximum TRT in each angular direction.

Rim radius: the lateral distance between the foveal center and the rim.

Maximum slope: the maximum derivative value in the region from the foveal center to the rim.

Parameters were estimated for all 24 angular directions and then averaged to obtain a single value per parameter and subject.

To locate the center of the fovea the methods described below were compared:

None: assume the center of the acquired scan as the foveal center.

Min: locate the foveal center at the A-Scan point of minimum TRT in the central 0.85 mm radius region.

Interpolation + min: resample the central part of the TRT map to a regular grid of 0.85 × 0.85 mm² and a 0.02 mm spacing using cubic interpolation. Then, locate the foveal center at the grid point with minimum TRT.

Smooth + min: resample the central part of the TRT map to a regular grid of 0.85 × 0.85 mm² and 0.02 mm spacing, and smooth it before locating the foveal center at the grid point with minimum TRT. We used the implementation of AURA Tools (foveaFinder.m function) [44] to smooth the resampled TRT map by applying a filter with a 0.05 mm radius circular kernel.

The main characteristics of the compared foveal pit mathematical models are shown in Table 2. Additionally, two smoothing methods were also applied with different degrees of roughness: moving average with five to 60 averaged samples, and local estimated scatterplot smoothing (LOESS) based on a second-degree polynomial with span in the range 1–50%. The smoothing was applied to each B-scan separately.

Characteristics of the compared mathematical models. The modelled region accounts for the part of the data that is modelled by each fit of the model. The number of parameters refers to the number of coefficients to be estimated in each fit.

^& The whole 2D total retina thickness (TRT) map is adjusted in one fit. ^$ The fovea is modelled radially using the foveal center as the reference. * The inner part of the B-scan (foveal center to rim) is fitted with two parameters while the outer part (the rim and beyond) is adjusted with three.

All the approaches were used to compute the parameters described in Section 2.3. However, the estimation of the CFT by the model proposed by Scheibe et al. [31] was equal to the raw estimation as the model uses the foveal center as a fixed reference point. Similarly, the model by Yadav et al. [28] uses both the foveal center and the foveal rim as fixed points and, therefore, did not affect the estimation of neither the CFT nor the rim height.