Automated Fovea Detection in Spectral Domain Optical Coherence Tomography Scans of Exudative Macular Disease

In macular spectral domain optical coherence tomography (SD-OCT) volumes, detection of the foveal center is required for accurate and reproducible follow-up studies, structure function correlation, and measurement grid positioning. However, disease can cause severe obscuring or deformation of the fovea, thus presenting a major challenge in automated detection. We propose a fully automated fovea detection algorithm to extract the fovea position in SD-OCT volumes of eyes with exudative maculopathy. The fovea is classified into 3 main appearances to both specify the detection algorithm used and reduce computational complexity. Based on foveal type classification, the fovea position is computed based on retinal nerve fiber layer thickness. Mean absolute distance between system and clinical expert annotated fovea positions from a dataset comprised of 240 SD-OCT volumes was 162.3 µm in cystoid macular edema and 262 µm in nAMD. The presented method has cross-vendor functionality, while demonstrating accurate and reliable performance close to typical expert interobserver agreement. The automatically detected fovea positions may be used as landmarks for intra- and cross-patient registration and to create a joint reference frame for extraction of spatiotemporal features in “big data.” Furthermore, reliable analyses of retinal thickness, as well as retinal structure function correlation, may be facilitated.


Introduction
Spectral domain optical coherence tomography (SD-OCT) is a noninvasive modality for acquiring high resolution, 3D cross-sectional images of the retina, and is today the most important ancillary test for the diagnosis and management of macular diseases [1]. Usually, serial OCT acquisitions are compared over time to determine disease progression and/or treatment response. However, as a result of acquisitions at multiple time-points, motion or imaging related registration errors commonly occur, that is, scans of the same eye at different time-points may be aligned incorrectly. Such registration artefacts may have a severe effect on the ability to perform accurate and reproducible analysis of subtle changes over time [2][3][4][5]. This problem can be overcome by computationally aligning sequential OCT scans using automated registration [2], or by labour-intensive manual alignment of the OCT scans or measurement grids.
In both cases, the fundamental requirement for registration of OCT volumes is the use of adequate anatomical landmarks. Further to the retinal vasculature, which has been used previously to register OCT volumes [2], the fovea centralis is a particularly important registration landmark. For example, correct identification of the foveal position is key for automated retinal thickness assessment using fovea-centered measurement grids such as the early treatment diabetic retinopathy study (ETDRS) grid and rotational alignment of the measurement grid in circumpapillary retinal nerve fiber layer measurements for glaucoma and multiple sclerosis. Due to its relevance for visual acuity, knowledge of the exact foveal position is also critical for studies of structure-function correlation [6]. However, given the immensely large scale of imaging data both in modern clinical practice and research, fully automated analysis methods are required for efficient OCT analysis. To our knowledge, computational detection of the fovea in SD-OCT has been limited to healthy or dry-AMD cases [7][8][9][10][11][12][13]. Thus, the major challenge of this work is to develop a detection method for the fovea in SD-OCT scans of exudative macular disease where the retina has been heavily deformed by the presence of fluid, severely altering the retinal anatomy.
In this article, we present a fully automated fovea detection method that is capable of accurately and reproducibly identifying the position of the fovea in cross-vendor longitudinal OCT scans of both normal and patients suffering from exudative macular disease, that is, cystoid macular edema secondary to retinal vein occlusion (RVO) and neovascular age-related macular degeneration (nAMD). In our method, we consider specific disease morphology to account for the differences between disease types. Our goal is to demonstrate the accuracy of the presented fully automated detection system and thus the systems' feasibility for detection of the foveal position in "big data."

Imaging Data.
For this study, 704 clinical SD-OCT imaging datasets from the Vienna Reading Center's (VRC) image database were used. 494 scans were selected from multicenter trials evaluating ranibizumab for central or branch RVO (clinicaltrials.gov identifiers, NCT01535261 and NCT01599650), and 210 scans from a multicenter trial evaluating ranibizumab for nAMD (clinicaltrials.gov identifier, NCT01780935). The study was conducted in compliance with the tenets set forth in the Declaration of Helsinki. The randomized clinical trials from which the scans were obtained were approved by the institutional review board of each participating center. All patients gave written consent for participation in the respective trial and all data was appropriately anonymized.
For method development, distinct training and testing image datasets were constructed. The training set, used to optimize the detection algorithms, consisted of 180 scans, divided into three equally sized groups representing branch RVO, central RVO, and nAMD, randomly selected from the baseline time-point where disease is most prevalent. The unseen testing set (used to validate the final algorithms) was comprised of 240 scans divided into the disease groups branch RVO (42 Heidelberg Spectralis, 33 Zeiss Cirrus, 5 Topcon), central RVO (53 Heidelberg Spectralis, 23 Zeiss Cirrus, 4 Topcon), and neovascular AMD (48 Heidelberg Spectralis, 32 Zeiss Cirrus), comprised of 80 scans each (all distinct from the training set), randomly selected to be inclusive of various acquisition time-points. In both the training and testing datasets, each scan was acquired from a distinct patient.
For all 420 scans, to provide an objective and standardized ground truth, the position of the fovea was manually annotated by expertly trained graders from the VRC using validated custom software [14]. The diagnostic criteria for the foveal center included (1) minimization of the retinal nerve fiber layer (RNFL) thickness; (2) presence of a foveal depression; (3) focal elongation of the photoreceptor outer segment signal, as described previously [15].
Inclusion criteria for dataset construction included acquisition from multiple devices (Zeiss Cirrus, Heidelberg Spectralis, and Topcon 3D OCT 2000) as well as baseline and nonbaseline time-points. Each OCT volume has average physical dimensions of 6 × 6 × 2 mm 3 ( , , ) with slice thickness ranging from 11 Figure 1 illustrates the three major stages comprising the International Journal of Biomedical Imaging automated fovea detection algorithm presented here, that is, image preprocessing, fovea type distinction, and fovea type based fovea position detection, which are discussed in detail below.

Fovea Detection Preprocessing. The flow diagram in
The preprocessing stage is shown in Figure 2. Firstly, motion correction in the plane is performed on the entire image volume using the local-symmetry estimation method described in [16] to remove motion artefacts caused by microsaccadic eye movement. This can be seen in Figure 2 in the left image showing the retina orientation flattened. Secondly, tilt correction, also described in [16], is performed in the B-scan plane reducing the horizontal tilt of the retina. The resulting motion corrected volume is denoted as V C .
Thirdly, denoising is performed on V C to reduce speckle noise using the block matching collaborative filtering approach described in [17], giving V NF .
On V NF , a kernel graph-cut segmentation algorithm [18] is applied to delineate the inner limiting membrane (ILM), as well as intraretinal cystoid fluid (IRF) regions, resulting in the segmented volume V G . Finally, computational complexity is reduced by masking V G . A 2D elliptical mask M, of size M = /2 and M = /2, is constructed in the XY plane centered at a statistically generalised mean fovea position, where and are and image dimension sizes. This center point has been computed by averaging manually annotated fovea positions from the training dataset, annotated by expert graders at the VRC, and was found to be within = 80 m and = 140 m of the scan center, based on their relative distances from the respective scan centers. As a result, the scan center is used, with all image information outside the mask area removed from the plane which is then propagated into 3D, resulting in a cylindrical region containing the fovea, V Gm . This cylindrical region is exemplified in the third image of Figure 2, showing the masked B-scan.
After this preprocessing stage, ILM delineation is performed again on V Gm to extract the cropped ILM surface for fovea appearance analysis.

Fovea Appearance Classification.
Given the distinct appearance of the fovea between healthy and diseased, as well as disease-specific variability, different methods for fovea position detection are required. In exudative maculopathy, pathologies that affect the fovea may be categorised into intraretinal cystoid fluid (IRF) resulting in foveal edema, IRF resulting asymmetrical foveal edema, subretinal fluid (SRF), and pigment epithelial detachment (PED). Thus, to further simplify the detection problem, we can generalise the fovea appearance into a small finite number of types based on their appearance when normal and pathological. Three main appearances of the fovea caused by the above described pathologic lesions have been reported previously [19] (Figure 3). In the first case ( Figure 3(a)), the fovea appears as a depression which we denote as a normal foveal depression (NFD). This is the appearance of the fovea when little to no disease is present. Cases with more severe pathology may be categorised into two major appearances, further exemplified in [20], where the fovea has been deformed by pathology such as IRF.
In minor foveal depression (MFD, Figure 3(b)), the fovea appears as a depression smaller than in a NFD case that has been raised in direction by macular edema. Finally, in absent foveal depression (AFD, Figure 3(c)), the fovea is not visible as a depression; instead, the ILM appears as a parabola, again elevated in the dimension by retinal edema below.
Automated Distinction of Foveal Appearance Types. Automated computational diagnosis of the NFD type examines the RNFL layer thickness using the clinical description of minimum RNFL thickness to describe the fovea in normal cases. This method has been previously described in detail in [19].
Distinction of the MFD type focusses on the unique morphology of this fovea appearance in the form of a minor depression that has been elevated in z-axis primarily due to IRF regions (Figure 3(b)). From the ILM segmentation in V Gm computed in the preprocessing stage (Figure 3(e)), maxima examination of this surface section is performed on a B-scan by B-scan basis. For the scan to feature a MFD, the masked surface must be comprised of 2D surface segments treated as a curve function featuring a global maximum and a further local maximum representing the two peaks seen in Figure 3(b). A peak is defined as a data point in the curve that is larger than its two neighbouring data points. Confirmation of maxima presence is performed by extension into the third dimension, where the required maxima configuration must be present across a contiguous distance of 150 m, given a mean fovea centralis diameter of 1.5 mm. A physical distance is used rather than number of B-scans due to crossvendor variable slice thickness. This is further exemplified in  Distinction of the AFD is based on the identification of a global maximum, across a contiguous 150 m distance, similar in appearance to a conical shape, as seen in Figure 3(f).

Fovea Position Detection.
From the fovea appearance classification stage, a given retinal SD-OCT volume is classified as featuring one of the three fovea appearances described in Section 2.3.1. In the fovea position detection stage, appearance specific detection functions have been developed to compute the fovea position, as described in this section.

NFD Fovea Detection.
In regard to anatomical features, we know that the fovea is the point at which the RNFL thickness is zero. Thus, in the NFD case, we can delineate and extract the two required surfaces (ILM and RNFL) using the graphcut based retinal surface segmentation algorithm described in [21]. In the NFD detection method, we are only interested in zero thickness in the fovea region, which has been masked and represented as V Gm . NFD fovea position detection is illustrated in Figure 4.
In the event that multiple zero thickness points are identified, the center of mass is computed and taken as the fovea position. This method is universally applicable to NFD scans.
International Journal of Biomedical Imaging  MFD Fovea Detection. The minimum thickness method for computing fovea position described in the NFD case can no longer be relied upon in the MFD case as this foveal appearance features a retina that has been deformed by pathology such as IRF. As a result, accurate RNFL segmentation is no longer reliable, specifically in the region of interest around the fovea where retinal edema is most prevalent. Thus, to locate the position where the RNFL is thinnest, we compute the distance between the ILM surface and the IRF causing the disruption in the retinal anatomy, as illustrated in Figure 5. Analysis of the manual fovea positions in the training dataset has shown that in all cases where the normal foveal depression was elevated by the presence of IRFs and creating the MFD, the B-scan on which the fovea annotation was performed corresponded to the B-scan where the distance between ILM and IRF boundary was thinnest in all corresponding cases.
However, due to the appearance and morphology of the deformed pathological fovea region, accurately computing the distance may be challenging. Thus, an additional preprocessing step is required, that is, ILM delineation of V NF in pathological scans.
Due to the graph-search based retinal surface segmentation algorithm described in [21] performing less adequately in pathological scans as opposed to normal scans, a proprietary ILM segmentation algorithm was developed based on the kernel graph-cut method described in [18]. In this multiregion approach, image partitioning is achieved using kernel mapping of the SD-OCT B-scan. Each B-scan is transformed implicitly by a kernel function in order to apply graph-cut formulation to the problem. In this case, we applied = 2 to describe the number of regions to segment, as well as a relaxed smoothing constraint to ensure that the ILM surface boundary is delineated as accurately as possible, as opposed to the smoothed surface from [21]. This ensures that the ILM surface used for RNFL thickness measurement in pathological cases is accurate. Vector ILM ( , ) describes this surface.
From V Gm , given delineated ILM from pathological scan and that the IRF below the fovea is labelled as the lowest intensity region within the retina delineated in V G , the distance between ILM and IRF is computed to calculate the fovea position.
Vector ( IRF ( , )) is computed from the segmented candidate IRF region from within V Gm comprised of the boundary points of the region. Furthermore, the minimum distance between ILM (upper surface, Figure 6) and IRF 6 International Journal of Biomedical Imaging  points (lower surface, Figure 6) is computed using pairwise Euclidean distance computation (1) between ILM ( , ) and IRF ( , ) (arrows, Figure 6): where is the vector of distances between the vector of points ILM (ILM surface) and IRF (IRF surface). The resulting pairwise distances are sorted in ascending order (ignoring the anatomically impossible zero thickness), choosing the shortest distance (thickest arrows, Figure 6) and the corresponding point from ILM ( , ). However, the volumetric characteristic is also important to consider as the retinal SD-OCT scan is a volume. Thus, similar to the appearance classification phase, the fovea detection algorithm is performed on every B-scan in the masked fovea region in the event multiple B-scans that have identical minimum distances. In this event, the center of mass is computed as the fovea position in the plane with position obtained from the ILM segmented.
AFD Fovea Detection. The major characteristic of the AFD is the parabolic appearance of the ILM, again due to deformation by IRF, and, as such, the same method used to compute the fovea position in MFD cases is applied. This is illustrated in Figure 6(a). Furthermore, as with the NFD case a similar correlation was found between the manually annotated fovea positions in AFD cases within the test dataset and distance between ILM and IRF boundary.

Statistical Analysis.
The performance of the developed algorithms was evaluated on the unseen validation dataset.
For categorical variables (i.e., fovea appearance type), the accuracy was descriptively analysed using confusion matrices and the area under the receiver operating characteristics (ROC) curve. Furthermore, the agreement between the automated and manual diagnosis was evaluated using Pearson correlation. For continuous variables (i.e., fovea position), the distance between manual and automated fovea positions was described as mean with 95% confidence intervals. Furthermore, the correspondence between manual and automated fovea positions was characterized using Pearson's correlation coefficient. The formal significance level was set at < 0.05.

Results
Implementation of the proposed method was carried out using MATLAB (Version R2012b, The Mathworks Inc.) on an Intel Core i7, 3.5 GHz, with 32 GB RAM.
Of the testing dataset described in Section 2.1, four central RVO (cRVO) and 6 branch RVO (bRVO) scans were excluded due to poor signal and image quality from acquisition. Thus, the validation dataset comprised 230 SD-OCT scans (76 cRVO, 74 bRVO, and 80 nAMD) SD-OCT scans. Table 1 presents the fovea appearance classification validation, comparing system results with ground truth fovea appearances labelled as either NFD, MFD, or AFD. Automated fovea appearance classification resulted in an 84%, 89%, and 88% correct appearance distinction for bRVO, cRVO, and nAMD, respectively, based on comparison with expert annotation. The area under the ROC curve (AUROC) was computed as 0.956, 0.949, and  0.938 for bRVO, cRVO, and nAMD, respectively. Furthermore, agreement between grader and automated appearance classification using the Cohen's kappa coefficient ( ) was = 0.741, = 0.804, and = 0.777 for bRVO, cRVO, and nAMD, respectively. Table 2

Discussion
In this article, we present a fully automated system for classification of the foveal shape and detection of the foveal position in SD-OCT volume scans with exudative maculopathy. Validation against manual ground truth provided by an experienced reading center demonstrated excellent accuracy of the automated system. Furthermore, our method showed applicability across nAMD and RVO diseases as well as across several prevalent SD-OCT devices.
Analysis of failure cases for bRVO and cRVO scans attribute incorrect appearance classification to incorrect delineation of the ILM surface, as a result of image/signal quality. As can be seen in Table 2, the largest proportion of error cases was MFD classified as AFD in both RVO types. Examination of these cases show poor signal and image quality in the B-scan plane within the retina as a result of acquisition/scanning artefacts and inability for either the layer segmentation described in [21] or our proprietary ILM segmentation to delineate an accurate ILM surface. As a result, the delineated ILM surface is not representative of the actual surface appearance, causing the incorrect classification of MFD as AFD. However, in nAMD cases, the majority of failures occurred when the system identified MFD as NFD. This is attributed to zero thickness computation by the system whereas a human grader has identified an elevation of the foveal depression as dictated by the MFD appearance guidelines.
In the second major contribution of this work, examination of the results of fovea position detection indicates a low absolute distance between the automatically computed fovea positions and manual ground truth for bRVO, cRVO, and nAMD, as well as on a per device basis. Overall, the most accurate automated fovea position detection was seen in bRVO, followed by cRVO, and then nAMD. Examination on a per device basis shows the fovea to be detected most accurately in Spectralis images for bRVO and Cirrus images in both cRVO and nAMD. Based on image quality, it would be expected that detecting the fovea within Spectralis images results in the highest accuracy; however, fovea appearance classification and position detection also rely heavily on the delineation of pathology, in this case, IRFs. As a result, the quality of the imaged pathology, which varies from acquisition to acquisition and is affected by other factors such as imaging artefacts and patient motion, will also play a role in the accuracy of the resulting detected fovea position. Comparisons of the automatically detected fovea positions against their manually annotated counterparts show high correlation (>0.9) for and fovea positions in all disease 8 International Journal of Biomedical Imaging types, except for fovea position in nAMD. The probable cause for this poorer correlation is likely due to the more variable disease features of AMD in comparison to RVO, thus resulting in more variable fovea position as detected by the automated method. Nevertheless, relative distances between the automatically detected fovea positions and their manually annotated ground truth counterparts remain low.
The examination of automatically detected fovea positions use manually annotated positions for evaluation; thus the accuracy of the ground truth must also be taken into consideration when evaluating the automated results. Not only is this task time consuming to perform precisely for human graders, particularly in a "big data" setting, but also the criteria used to determine the position will also be affected by human subjectivity. Previous studies revealed an mean interobserver variability of 71.63 m for the foveal position but also a mean distance between the true and device center point of 290.9 m [15]. Thus, the agreement between human experts is higher compared to our automated system; however, our system accuracy is greater than that obtained by the devices specified in [15], at a mean of 195.5 m, illustrating the applicability of our method in comparison to other center point detection algorithms and as a more practical alternative to manual center point plotting in a "big data" environment. Furthermore, B-scan spacing must also be considered as the disparity in the number of B-scans from Spectralis, to Cirrus and Topcon varies up to a factor of 1 : 5. As a result, computation of distance between system and ground truth will be affected. For example, a fovea position misaligned in dimension by a single B-scan in Cirrus (200 B-scans) would result in an error of ∼30 m, whereas the same misalignment in Spectralis would result in an error of ∼125 m. Thus, such disparity will affect human grader ground truth, and, by extension, system result validation, and given that the test dataset used here features Heidelberg Spectralis scans in the majority (∼60%, 143 of 240), explaining the higher distance error reported by the system presented here.
To the best of our knowledge, this is the first fully automated method to locate the position of the fovea directly within retinal SD-OCT scans, independent of whether they are highly pathological due to the presence of cystoid macular edema secondary to central and branch RVO, and neovascular AMD, or within the healthy retina [7] or only via a fundus photograph [8][9][10][11][12][13]. Future work would concentrate on extension of the IRF segmentation to incorporate a machine learning approach that distinguishes IRF from SRF, allowing for targeted IRF delineation as a feature for computing the fovea position. In addition, the relationship between the fovea region and the thickness of other retinal layers may be explored and used as further anatomical features for detection.
In summary, we have presented a fully automated approach to detect the fovea within healthy and diseased SD-OCT scans of the macula. This enables the use of the fovea as a key landmark in the construction of a population reference frame to identify and extract key spatiotemporal features from a large patient population comprised of different time-points, devices, and imaging modalities. Furthermore, being the functional center of vision, the fovea is crucial for performing analyses of retinal structure/function correlation [22,23].